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MODELING OF THIN-FILM GaAs GROWTH 


By 

John H. Heinbockel* 

SUMMARY 

A potential scaling Monte Carlo model of crystal growth is developed. 
The model is a modification of the solid-on-solid method for studying crys- 
tal growth in that potentials at surface sites are continuously updated on a 
time scale reflecting the surface events of migration, incorporation and 
evaporation. The model allows for B on A type of crystal growth and lattice 
disregistry by the assignment of potential values at various surface sites. 
The surface adatoms are periodically assigned a random energy from a 
Boltzmann distribution and this energy determines whether the adatoms evapo- 
rate, migrate or remain stationary during the sampling interval. For each 
addition or migration of an adatom, the surface potentials are adjusted to 
reflect the adsorption, migration or desorption potential changes. 

INTRODUCTION 

Numerous methods have been applied to obtaining thin-film, single crys- 
tals of GaAs, including free-standing wafers, peal films removed from a 
single crystal substrate, and films grown on lightweight substrates. The 
most promising method is a version of the latter technique called 
"graphoepitaxy. " It is generally known that overlayers of crystalline mate- 
rials deposited upon smooth microcrystalline substrates tend to be more or 
less randomly polycrystalline. The absence of long-range order in the 
microcrystalline substrate is reflected in the absence of long-range order 
in the over layer. The basic concept of graphoepitaxy is that, by introduc- 
ing an artificial surface relief structure having a long-range order on a 
microcrystalline substrate, long-range order can be induced in an overlayer. 
In other words, a crystalline film can be grown on a raicrocrystalline sub- 
strate . 

*Professor, Department of Mathematical Sciences, Old Dominion University, 
Norfolk, Virginia 23508. 


The degree of crystalline order achieved during a growth process will 
be controlled by the adsorption, nucleation and lateral growth behavior in 
the first few deposited layers. The parameters which will affect the crys- 
tal growth are: The deposition rate, the surface temperature, surface dif- 

fusion, surface defect density and lattice registry of the system. We 
present a Monte Carlo solid-on-solid (SOS) computer simulation which uti- 
lizes a potential scaling technique (ref. 1) over a 20 x 20 array of sites. 

Although numerous Monte Carlo models for crystal growth exist (refs. 2-12), 
the approach developed herein is a more physical model in that the events 
which occur at each site are constrained by the surrounding potential field 
and the thermal energy fluctations associated with a given substrate temper- 
ature. Also, the method of ordered statistics is utilized to construct a 
time scale of events compatible with computer times in order that the simul- 
taneous changing and updating of site potentials can be done in a reasonable 
amount of computer time and still allow the model to simulate thin film 
growth in a physically realistic manner. 
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Boltzmann constant (8.62 x 10“ 5 ) eV/*K 
absolute temperature (°K) 
scale factors for crystal orientation 
crystal orientation factors 
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deposition rate (nm/sec) 
evaporation rate 
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random energy (eV) 
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sample size 

site numbers 
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cumulative distribution 

distribution of ordered statistic E(n) 
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crystal orientations 

heat of sublimation 

heat of adsorption for single adatom 

diffusion activation energy (eV) 

lifetime (sec) 

mean time between hops 

output snapshot or stop action time interval (sec) 
incorporation energy 
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kink site potentials 


Mie potential parameters 
diffusion coefficients into and from bulk 
sticking coefficient 

period of vibration of surface atom (sec) 
parameter for surface diffusion 

MODELING OF CRYSTAL GROWTH 

The SOS Monte Carlo Model 

In this model we consider a simple cubic (SC), body centered cubic 
(BCC) and face centered cubic (FCC) crystal lattice. It is required that 
each occupied site be directly above another occupied site and so the name 
solid— on— solid (SOS) model. This model is characterized by an array of 
interacting columns of varying integer heights with respect to some orienta- 
tion such as the (100) (111) or (110) crystal planeo. The terrace-ledge- 
kink (Kossel) Model (refs. 13-14) is illustrated in figure 1. The model 
employs a 20 x 20 square array upon which columns are constructed. Adatoms 
are deposited upon the surface in a random fashion where they are free to 
migrate, remain localized, diffuse into the bulk (incorporation) or diffuse 
from the bulk to the surface, or evaporate. These changes alter the 
stacking heights of each column as well as producing new potentials at and 
in the neighborhood of the surface sites involved in the process of surface 
adatom interaction. 

Each surface site is located at the top of the stack of adatoms from an 
arbitrary row i and column j of a 20 x 20 array. The physical 
constraints which determine the temporal behavior of every adatom located at 
the surface of each column from an arbitrary site (i,j) is based upon the 
interaction potential that the surface adatom has with its nearest neighbors 
and the rest of the solid. Spatial disregistry, that normally occurs due to 
size differences between adsorbate and adsorbent atoms, is accounted for by 
changes in the interaction potential. In figure 2(a), the interaction 


(1) „ (2) 
ks ’ ks 


* * 
m , n , d> , R 
o’ o * r ’ 


D-, D + 
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potential across a perfectly homogeneous surface is depicted as uniformly 
changing from site to site. Figure 2(b) illustrates a typical interaction 
potential across a heterogeneous surface. We developed a set of "rules" 
whereby the columns of the SOS model interact by assigning values to the 
potential energy changes associated with th% processes of adsorption, 
migration and desorption. 

Potential Scaling of Adatoms 

The rules, by which the columns of the SOS model interacted, were gov- 
erned by the following ideas relating to the potential energy and potential 
energy changes associated with the adsorption, migration, or desorption of 
adatoms from an arbitrary row i and column j of an 20 x 20 array. Ener- 
gies associated with an arbitrary site (i,j) were defined as follows: U Q 
“ U Q (i,j) — the potential energy at a site because of surface bonding and 
crystal structure; <|> 0 — the potential energy change at site (i,j) because 
of the deposition of an adatom (assumed the same for all sites); -W£(i ■ 

1, ..., 8) — the potential energy changes at neighboring sites when an adatom 
is deposited at site (i,j); E(i,j) — tha random surface energy associated 
with site (i,j) and time interval At s ; U(i,j) = U Q (i,j) + E(i,j) — the 
total energy associated with site (i,j) during the time interval At 8 ; 

U e — the evaporation potential; and U m — the migration potential. All of 
the above energies are measured in electron volts. 

We developed a Monte Carlo computer simulation of crystal growth by 
developing rules that determined the SOS kinetics of condensation, evapora- 
tion or surface migration of adatoms. These rules led to a consistent and 
physically reasonable description of the fundamentals associated with crys- 
tal growth. We first considered the adsorption of a thermally accommodated 
adatom onto the surface at some general site where the potential at this 
site was changed and, simultaneously, potential energy changes at all of the 
neighboring sites occurred. In Table 1, the potential energy changes are 
depicted by the mnemonic mask. The center of this mask is placed over the 
site (i,j) to illustrate the changes to be made in the potential at the 
central site as well as the potential changes in the surrounding neighboring 
s ites . 
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The potential changes in the case of desorption of an adatom from the 
central site are again depicted with the mask of Table 1, with the opposite 
signs on the potential changes. The case of surface migration was treated 
as a desorption from a site (i,j), followed by an adsorption at a nearest or 
second nearest neighbor location, together with the correct potential mask 
changes associated with each process. The neighboring migration site was 
determined by a weighted random walk to one of the unoccupied neighbor 
sites . 

Table 1. Potential energy changes associated with central site (i,j) 
and neighbor sites due to deposition of an adatom (for a 
100 orientation). 


-W7 ■ -w7(i-l, j-1) 
-we " -W6(i,j-1) 
-W5 ” -w5(i+l,J“l) 


-w8 “ -w8(i-l,j) 
<J>0 * <}> 0(i, j) 

-W4 ® -W4(i+l,j) 


-Wi * -wl(i-i, j+1) 
-W2 ** -W2(i,j+1) 
~W3 » -w3(i+l,j+l) 


The Monte Carlo simulation of crystal growth involved a random deposi- 
tion of thermally accommodated surface adatoms during a time interval At. 
These deposited adatoms changed the potential energies at the random surface 
sites under consideration. The values assigned to the central potential 
change «#> 0 and neighboring potential changes -W£, i * 1, ..., 8 dic- 
tated the new potential energy values when an adatom was deposited or re- 
moved from a site. In this way each surface site had an energy barrier to 
translation or evaporation, represented by a potential well. 

Figure 3 illustrates the potential changes that occur along a lineal 
section of a homogeneous surface upon the adsorption of a single adatom at 
site (i,j). Note that the potential increases by an incremental amount 
<J> 0 , and the adjacent sites decrease in potential by an incremental amount 
^ . This represents the actual physical condition that the adsorbed adatom 
requires less energy to desorb or to migrate as compared to the original 
surface adatom which was surrounded by all its nearest neighbors. Note also 
that the deeper potentials at the neighboring sites reflect the increased 
energy necessary to desorb an atom from these sites due to the increased 
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coordination or ligancy created by the adatoms. Figure 4 illustrates the 
potential variation at ati arbitrary site (i,j) as nearest neighbors are 
progressively added onto a ( 100 ) surface orientation. 

The term epitaxy means "an arrangement on," and is used to denote the 
growth of one substance upon the crystal surface of a foreign substance. 

The term autoepitaxy refers to the oriented growth of a substance ~nto it- 
self, and hetroepitaxy is the growth over another material. Autoepitaxy 
requires that the initial potential U 0 (before adsorption) be recovered 
when the adatom has all its neighboring adatoms surrounding the central 
site. 

In assuming values to the potential changes ^ and i“l,.,. 8 , we 

must take into account the type of crystal structure and orientation we are 
trying to simulate with our SOS model. Consider figure 5 which illustrates 
the GaAs fee structure. For growth on the (100) face we can set up a cor- 
respondence between a central site, the nearest neighbor sites, second near- 
est neighbor sites, and the adatom potential changes for the mask in Table 1 
(i.e., “ <{> 2 , “ <}>i, etc.). Similarly, we can sat up the correspondence 

illustrated in figures 5(b and c) for the (111) and (110) orientations and 
we can construct an appropriate mnemonic mask. 

In figure 5 we must choose 9o> < t , i> < f , 2> < f > 3 Buch a wa y that > when ttie 
first level of adatoms covers the surface, then the potential distribution 
must return to its original value. We will require that adjustments be made 
in the potential energy changes during the transition from heteroepitaxy to 
autoepitaxy. Here we let a negative sign denote an attractive potential. 

By simply adding adatoms to a surface it is readily verified that the poten- 
tial changes must adhere to the rules given in Table 2 if after one layer 
the potential energy returns to its initial value. 

For heteroepitaxy we require that an adjustment be made in the central 
site potential change to reflect the potential energy differences of the 
materials involved. The potential energy changes <p 0 , , <p 2 > $3 can be 

different for the substrate and growing material. For the substrate materi- 
al we could use the depth of the surface potentials and migration levels to 
stimulate a variety of surface morphologies. In this model we envision a 
flat substrate as a periodic lattice structure where each lattice site is 
a potential well. The substrate can vary from flat to rough and the 
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potentials adjusted to reflect various surface preparations. For an ideal 
ly flat substrate we assume that the depths of the potential wells are 
uniform, given by U . After one layer of growing material covers the 
surface, the potentials at each site are assumed to convert to the auto- 
epitaxy potentials U 0 . In order to make this transition we assume that 
8 

$ * I to. + (U - U )T. . where T. . is zero if the height h.. at 

position (i,j) is greater than or equal to one, and is one in the case 

where h.j “ 0. Thus, if an adatom is deposited at a first layer site 

(i,j), we adjust the potentials at this site by the relation U - U , in 

o 08 

addition to the mask potential changes at the surface sites as this produces 
the desired change that hetroepitaxy produces in the value of the poten- 
tial. 

The energy behavior at each surface site is monitored over a sampling 
interval At s , every sampling interval. The simulation of thermal energy 
fluctuations is done by random number generation. For each site (i,j) we 
generate a random energy E(R) and add it to the interaction potential 
U 0 (’ j,j) ?aa obtain a total energy of 

U(i, j) = U Q (i,j) + E(R). (1) 

This eneigy is then compared to energy barriers for desorption, surface 
migration and incorporation. If the total energy exceeds one of these bar- 
riers, we allow the adatom to proceed accordingly. If the total energy is 
less than the lowest barrier, then the adatom remains localized. In sum- 
mary, our Monte Carlo procedure entails the generation of a random energy 
for each site and calculating a total energy U, if this value U is such 
that: 

(a) U < U , the adatom remains localized; 

(b) U < U < U., surface migration is allowed to occur; 

(c) U. < U < U , surface migration or incorporation into the 

bulk is allowed to occur; 

(d) < U, evaporation occurs. 


8 


During each time interval At s , a random energy E(R) was assigned 
to each of the surface adatoms. We let 

U(i,j) “ u 0 <i,j) + E(R) (2) 

denote the total energy possessed by a surface adatom at a site (i»j) during 
this time interval. This total energy is the sum of the potential energy 
U 0 due to the lattice structure and a random energy E from a modified 
Boltzmann distribution to be discussed in the next section which character-” 
izes the random surface energy. When U was less than some material-de- 
pendent migration level U m , the adatom remained stationary at the 8ur‘» 
face site. If U m < V < U e , surface migration by random walk was allow- 
ed to occur. If U was greater than the evaporation potential U e , the 
adatom was removed from the site and for < U < U e incorporation or 
migration was assumed to occur. 

The rate of impingement of adatoms upon the surface was independent of 

the surface configuration. The rates associated with the evaporation and 

migration of adatoms depended upon the potential barriers U e and U m 

and also upon the values assigned to the potential changes and 

— W£ , (i - 1, ...» 8). Th^se later potential changes had to take into 

account the type of crystal structure and orientation of the growth we were 

trying to simulate with the SOS model. In figure 5(a), for growth on the 

(100) face, we sec up a correspondence between the central site, the nearest 

neighbor potentials <f) ^ , second nearest neighbor potentials <j> 2 , and the 

adatom potential changes for the mask in Table 1 (e.g., w 1 * $2’ w 2 * <J>i)* 

In a similar manner we were able to set up the correspondences illustrated 

in figure 5(b) and (c) for the (111) and (110) orientations. In Table 2, we 

selected the relation between the neighbor potentials 4> 0 , , <j> 2 , <j> 3 in 

such a way that when the first level of adatoms covered the surface, the 

potential distribution returned to its original value. This produced the 

constraint conditions on the neighboring potentials which are illustrated in 

Table 2. We assumed that # - AHads and were left with decisions on how 

o 

to assign the , (J , 2 ,<j>3 values. For the (100) orientation we let $3 “ 0 
and were left with having to assign values to ® ne P 0SS iEle choice 

was to assign <J>^ a value based upon nearest neighbor bond strength and 

then calculate the <J>2 value based upon the constraint. 
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Alternatively, we let ^ denote the change in the nearest neighbor 
potentials due to the addition of an adatom to the surface and let ( t > 2 ,< J ) 3 
denote the second and third nearest neighbor potential changes. We assumed 
that <f> 2 “ a 2^1 and 4*3 " a 3^1 w here a 2» a 3 are scale factors which are less 
than one. This allowed us to define the crystal orientation factor 5 as 


2 + 2c»2 

3 + 3a 2 

1 + a 2 + 2a 3 , (110) 


, (100) 
, (111) 


(3) 


which takes into account the different crystal orientations. We also defined 


( 1 ) 

ks 


( 2 ) 

and after and the capture of an 


the kink site potentials before 

adatom as ■ U - ■ U *• £<f>i. (Note that $ ■ 2C4>i.) Note 

ks o 1 ks o TA o 

that the values assigned to the mask potential changes are not necessarily 

the same for the different orientations: for example, the <j>i , <f> 2 > 4*3 values 

for each case in Table 2 could have different values. 


We can assign arbitrary values to the neighbor potential changes 

> ^2 >$3 as l° n S as these values satisfy the constraint that <f) 0 =* 25<j> 1 . 

If we arbitrarily assign values to a 2 » a 3 than we can solve for ^ and 

consequently <f> 2 , and <J>g. Instead of arbitrarily assigning values to a 2 

and cu, we examine a (m , n ) Mie potential curve with m < n (ref. 15) 

3 o o o o 



Here 




the distance at which <J> obtains its minimum value of 
If we examine the potential values at various distances 


R 


vr 

a„ , — a, 
2 


V3” V5 

— a_ , — a_ , with 
,2 ° 2 ° 




then the ratio of the Mie poten- 


tial values at these distances can be used to approximate the values of a 2 
^ 2/^1 and a 3 " ^ 3/^1 f° r different crystal orientations and different m Q , 


n n values. 
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Table 2. Potential changes for addition of an adatom to an arbitrary site. 


Crystal 

Face 


( 100 ) 


( 111 ) 


( 110 ) 


Potential Changes 


Relation Between 
Neighbor Potentials 

for Addition to 
Arbitrary Site 



-0)? 

-*8 

"“I 



““6 

^ 0 

-u 2 



-0) 5 

-«4 

- w 3 



CM 

-e- 

1 

1 

-♦l 

"*2 

<|>Q a 4$! + 4<|>2 


-♦l 

“$0 

-h 



-*2 

~*1 

— 1 

CM 

-e- 

1 


- 

~$2 




-♦2 

-<h 

-<*>1 

“+2 

4> o " 6< h + 6< ^2 


“*1 

*0 

“♦l 



"<j> 2 

-+l 

“<h 








“^3 


“*3 

<f> 0 a * 24 > 2 + 4<J»3 

“*1 

h 

“<<>1 



-*3 

-<t> 2 

-<j> 3 


Distances to 
Neighboring Sites 




Energy Distribution and Time 


During each sample interval At s> an atom in an arbitrary site (i,j) 
has the total energy 

U(i,j) “ U o (i,j) + E(R) 

where E(R) is determined by random number R from the Boltzmann distribu' 
tion 

f(E) » X exp [-XE], E > 0 (5) 

where X ** 1/KT is the parameter of this exponential distribution. The mean 
energy of this distribution is 

“ 1 

<E> = / E f ( E ) de - KT * - (6) 

o X 

and the cumulative energy distribution is given by 
E 

F(E) =* / f(E)dE = 1 - exp [-XE] . (7) 

o 

A random variate E can be generated from this distribution by using the 
inverse function associated with the cumulative distribution. For R a 
uniform random number between 0 and 1 and with R *= F(E), the inverse func- 
tion gives the random energy 

S = E(R) = - KT in (1 - R) (8) 


so that (1) becomes 

U(i , j ) - U (i,j) - KT In (1 - R) (9) 

o 

The residence or stay-time of an adatom on a surface is given by the 
Frenkel equation (ref. 16) 

r = x exp [-X AHads] (10) 
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where AHads is the heat of adsorption and T 0 is the period of vibra- 
tion for the surface adatoms (t 0 ~ 10“ 12 sec). Ideally then, the most 
physically real sampling time corresponding to changes in vibrational energy 
and, therefore, changes in U(i,j) is to choose At * T . Computer costs 
and time, of course, prohibit the extensive amount of computations that 
would be necessary to sample 400 sites 10 12 times each second. In order to 
circumvent this difficulty, the method of ordered statistics is applied. 
Essentially, most of the time-dependent energy variation at a particular 
site results in insufficient thermal energy for adatom movement and adatoms 
remain localized over most of the sampling interval. Since this large time 
of atomic localization is not important to the actual thin film growth, only 
the fraction of the sampling interval that movement does occur need be con- 
sidered. Thus we desire that fraction of the time that the site energy is 
in excess of the minimum activation barrier for adatom activity. 

Let Ep E 2 ,...E n denote n random samples from the exponential 
distribution (3) and let E^ lN , E(2)'**** denote the ordered arrange- 

ment (from low to high) of the n random samples with E^_^j < E^^ for all 
i=2,3, . . , ,n. The probability distribution of the largest ordered statistic 
E^ * max {Ej, E 2 ,..., E^} is given by (ref. 17) 

g(E) = n[F(E)] n_1 f(E), 0 < E < « (11) 

where f(E) and F(E) are given by equations (3) and (5). The cumulative 
frequency distribution is given by 

E 

G(E) = / g(E)dE « (l - e" AE ) n (12) 

o 

To generate a random variable E^ from this distribution, a uniform 
random number R, with 0 < R < 1, is generated such that 
the inverse function gives the random energy 

E, v « - KT An (1 - R 1//n ) (13) 

(n) 
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which can be compared with equation (6) . Note that for large values of n, 
we can approximate the random energy E( n ) by 

E, . « - KT *n (- — *n r). (14) 

(n) „ 

If for example, T ■ 500 K and the activation energy for diffusion is 
Q d ■ 0.7 eV, then the mean time between hops is x^ * 10“ 12 exp [XQ d ] ** 10” 5 
sec or, in a sampling interval of At s = 10” 5 sec, a single hop occurs. 

Any smaller sampling interval is not necessary because no movement occurs. 
Any larger sampling interval would result in multiple events for a single, 
adatom and the less physical the model becomes. Figure 6 illustrates the 
probability distribution of the ordered statistic with n ■ 10 7 , 10 8 , 

10 9 at T - 500 K. 

The minimum activation energy for diffusion (ref. 18) determines the 

sampling interval and therefore the number of random samples, n. We make 

the following assumptions concerning the activation barriers for adatom 

activity (see figure 3): U = AHsub, U = 0, <J> = AHsub - AHads, U = 

J o e o m 

-AHads + , x^ » T q exp (XQ^), then the mean number of hops in At g sec is 

given by At g /x an d n = At g /X Q . Let At Q denote the output "snapshot" time 

interval where we perform a stop action and view the surface. The number of 

surface scans during this time interval is given by At /At . For example, 

o s 

if At =0.1 and At = 10“** we would scan each of the 400 surface sites 10 3 
o s 

times and generate 4(10 5 ) random energies E^ from the modified Boltzmann 
ordered statistic probability distribution. 

Random Walk and Incorporation 

For a fixed uniform deposition rate R^j, adatoms are deposited at 
random positions on the surface based upon the value assigned to At s . 

For small At s , we must wait for some multiple of this time before adding 
a single adatom. Each time an adatom is deposited, or removed from a site, 
the potentials surrounding the site are updated. In the case of surface 
migration, an adatom has sufficient energy to migrate and we treat migration 
as an evaporation followed by a deposition at a neighboring site. The 
availability of more than one site for migration is another decision which 
is made according to the ligancy or coordination number associated with the 
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available sites. From an interaction potential perspective, this is physi- 
cally reasonable since the site that has a deeper well (more attractive) 
will have an energy barrier to migration that is smaller, thus having a 
higher probability for migration to that site. The experimental evidence to 
support this assumed behavior is sizeable (ref, 19). Each unoccupied near- 
est neighbor and second neighbor site is given a weight which is the ligancy 
if an adatom had random walked to that site. These weights are then normal- 
ized and a weighted random walk to one of these sites is determined by a 
random number. If all nearest neighbor and second nearest neighbor sites 
are occupied, the adatom is assumed to jump up onto the next level at an 
unoccupied nearest neighbor site with equal probability. 

In the case IL < U < U e> an adatom is considered to be either incorpo- 
rated or to surface migrate by weighting the two possibilities according to 
their relative probabilities. In the case of incorporation an adatom is 
removed from the surface in the same way an adatom desorbs. It is assumed 
that the bulk vacancy concentration is sufficient to receive the adatoms and 
therefore the adatom just disappears from the surface site. The excess 
energy after an event is assumed to be reabsorbed into the thermal energy of 
the solid. 


Generalizations 

Various modifications and extensions of the ideas presented in the 
previous sections will make the model more general. Some of the modifica- 
tions will be presented in this section as these reflect modifications of 
the SOS model. Three-dimensional model concepts which differ from the SOS 
model concepts will be presented in a later section. 

As adatoms are deposited upon the surface we will label them as "A- 
type" or "B-type" where the A-types represent substrate adatoms and the B- 
types represent the growing material. As the vertical growth increases the 
fraction of A's mixed with the B*s is allowed to decrease. In this way we 
can simulate adatom diffusion through the growing film to the surface by 
randomly depositing substrate adatoms on the surface, at a rate controlled 
by the diffusions equations for the adsorbate/adsorbent system. Also, 


instead of one potential change mask we will inti/ouuce four such masks to 
reflect the potential changes associated with the v, position process of: A 

on A, A on B, B on A, or B on B type of interactions. This additional 
complexity increases the bookkeeping to keep track of where the A and B^ type 
adatoms are located. 

By judicially choosing the A on A, A on B, B on A, and B on B potential 
changes we can use the SOS model to simulate A-B type crystal growth. By 
assigning an initial substrate of an A-B checkerboard pattern we can make 
the A on A interaction subject to a high probability of migration of A 
adatoms to a nearest neighbor B-site. Similar considerations hold for B on 
B interactions. The nature of interatomic potentials (ref. 20) can also be 
varied from site to site. 

Description of the Computer Program 

A flowchart of the computer program is given in figure 7 and the 
FORTRAN computer program is presented in Appendix A. The program is lib- 
erally spiced with comment statements to help the reader. An attempt was 
made to make the program modular in character in the event extensive revi- 
sion was needed. The following is a brief description of the subroutines: 

Program Crystal 

Here parameters are read in and other variables are initialized and, 
before the program actually runs, a printout of all initial values and para- 
meters is made. 

Subrouti n e SETUP 

Assumes a (m Q ,n ) Mie potential and calculates <j>^, a 2 , 013 given the 
initial value for <j> 0 . 

Subroutine FACE 

Calculates the potential changes associated with a (100), (111), or 
( 110 ) crystal orientation and distances to nearest neighbor site. 

Subroutine INITIAL 

Initializes and prints out the substrate geometry and assigned values 
of the potentials at each site. 
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Subroutine DXFF 


A weighted random walk surface diffusion i<* simulated by an evaporation 
followed by a deposition at the nearest neighbor site. 

Subroutine PBOND 

Calculates the weights associated with the random walk processes on one 
of the crystal orientations. 

Subroutine ADATOM 

Updates the potential sites at i,j and the surrounding neighbor sites 
when an adatom is added to site (i,j). 

Subroutine SUBATOM 

Updates all neighboring potentials when an adatom is removed from a 
site (i,j). 

Subroutines EDGE, EDD, XYX, CORRECT 

The geometry of the 20 x 20 square array assumes periodic boundary 
conditions as illustrated in figure 8. This 20 x 20 array is embedded into 
a 24 x 24 square array and addition or removal of adatoms along an edge, or 
migration across an edge, requires that the outer border of the 24 x 24 
array be updated to reflect the periodic conditions. These periodic condi- 
tions are maintained by the subroutines EDGE, EDD, XYX. and CORRECT. 

Subroutine STOFIN 

Starts and finishes a computer run. This subroutine analyzes the de- 
position rate and scan time and deposits adatoms on the surface, if re- 
quired. Surface scanning of each site is performed and computer output is 
printed. 

Subroutine URNS 

Calculates random sites and type of adatom (A or B) to be deposited on 
the surface. 

Subroutine UPDATE 

Scans the top layer of the surface and counts the size and frequency of 
clusters. Also calculated are coverage for the two highest layers, rough- 
ness factor, and number of A and B type adatoms on the surface. 
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Subroutine MIXUP 


Takes the vector (1,2,3, ... ,N) and randomizes the position of the 
integers to form a vector (u^ ,U 2 >u 3 , . , , ,u^j) which is some rearrangement of 
the initial vector. 

Subroutine GROW 

Optional computer graphics of output data. 

COMPUTER RESULTS AND DISCUSSIONS 

Various computer experiments were performed with the model. These 
experiments are listed in Table 3. The computer experiments were performed 
for a (100) fee surface to assess the physical behavior of the model. 

Because there exists no real data on semiconductor crystals for the 
interaction and potential energy changes associated with adatom movement, we 
chose the following set of nominal parameters: 

AHads ** 1.7 eV 

AH sub - 3.87 eV 

D_ z (Ge/Ge) =*7.8 exp (-2.98 X) cm 2 s“* 

D_ z (Si/Fe) ® 0.44 exp (-2.09 X) cm 2 s“* 

D+ z (Fe/Ge) “ 0.13 exp (1.08 X) cm 2 s" 1 

Q d = 0.7 eV 

S = S£ “ 1 

which represents the best available data for a Germanium type system. 

In experiment 1, we examined the surface migration of a single adatom 
which performed a random walk to the nearest neighbor sites whenever the 
random energy associated with a single scan dictated such behavior. The 
number of migrations as a function of time is linear and varies exponential- 
ly with temperature. Figure 9 shows a plot of diffusion coefficient vs. 
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Table 3. Computer experiments 



Experiment 

Deposition Rate 
(cm’^s"^ x 10*5) 

Temperature 

(K) 

Run Time 

1 . 

Surface diffusion of a 
single adatom 

0 

400, 500, 600 

2 

2( 

a) Nine adatoms in a row 

0 

400, 500, 600 

2 

(b) Nine adatoms randomly 
distributed 

3. Thin film growth with 

0 

400, 500, 600 

2 


defect site 

0.5 

400, 500, 600 

2 

4. 

Thin film growth with 

0.25, 0.5, 1, 
[2.5 (10-3)] 

500 

2 


variable deposition 

t>ULv< 

400 

40 

5. 

Thin film growth with 
variable substrate 
temperature 

0.5, 1.0 

300, 400, 500 
600, 700 

2 

6. 

Annealing of thin film 
growth 

0 

600, 700 

6 
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inverse temperature as calculated from the computer model. The linear be- 
havior of the diffusion constant follows from the equation 

j i a 2 

D - — - exp (-AQ ) (15) 

T 

O 

where a 0 is the jump distance and A ■ 1/4 for a (100) face. The slope 
of the Arrhenius plot yields the activation barrier Q<j " 0.7 eV which is 
our input condition. This provides a self consistent check on the physics 
of the computer model. In figure 9, the number by the circles is the 
average number of migrations for each temperature given. The results in 
experiment 1 are for a flat surface. In the case of nonuniform surface the 
mean diffusion length and migration frequency would be substantially reduced 
due to the lower probability of escape from kink sites, steps and other 
defects . 

In experiments 2(a,b) we examined the clustering of lineal and randomly 
dispersed adatoms as a function of temperature. In both cases the adatoms 
tended to seek out the most stable configuration in that each adatom ulti- 
mately tried to maximize its number of nearest neighbors and hence form a 3 
x 3 array. In the first case 2(a) the lineal adatoms were all connected and 
the probability of adatoms breaking away from the row increased with temper- 
ature. The adatoms tended to form stable clusters with the migration of 
single adatoms along & step being the dominant form of motion. For this 
experiment a 3 x 3 array was the most stable cluster. In the second case 
2(b) the dispersed adatoms performed random walks and collided to form 
dimers, trimers, and eventually a nine adatom cluster. Figure 10 is a 
graphic display of the sequence of events as a function of time. At some of 
the temperatures the 3x3 array was not completely cchieved. However, 
given sufficient time, these configurations eventually became a 3 x 3 array. 

The question of cluster growth is perceived to occur by either adatom 
capture or by actual enmasse motion of one of the smaller clusters into a 
larger cluster. The subsequent collision and reorientation of the adatoms 
to registry with the second cluster was assumed to occur. In the formula- 
tion of this model no consideration was given to the motion of whole 
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clusters, but as these computer experiments show whole clusters can move by 
individual adatom motion at the periphery of the cluster resulting in a net 
motion of its center of mass [see figures 10(b,c)]. 

In experiment 3, a generalized point defect was modeled by adjusting 
the potential to be very large at a single point and its four nearest neigh- 
bors sites. After one layer covered this site only the central site was 
allowed to have a larger potential and after a second layer the potential 
was allowed to return to the normal value of that of its neighbors. Figure 
11 shows the prescribed potential changes for the first two layers of 
growth. After the secorl layer covered the trap, the potential was allowed 
to return to its normal value. Figure 12 illustrates the effect of a trap 
and the resulting growth around the trap for a constant deposition rate and 
several surface temperatures. Many interesting phenomena occurred in our 
study of a trap and we will pursue this case in more detail in a later re- 
port. For the time being, a brief description of the observed phenomena 
will have to suffice. At low temperatures the surface adatoms initially 
adsorbed do not statistically interact with the site and ordinary homogene- 
ous nucleation and growth occurs uniformly over the surface. As adatoms 
impinge upon the trap site heterogeneous nucleation occurs and growth is 
much more rapid in this vicinity. Further the vertical growth in the vicin- 
ity of the trap is initially larger because the potential at the trap site 
is still lower after being covered by the first layer which is an island 
upon which impinging adatoms can migrate and find this lower potential. It 
is thus likely that the growth in the vicinity of the tsf&p would be 
dominated by this effect. At higher temperatures the major growth occurs by 
way of vapor phase transport as opposed to surface defect density (ref. 21) 
and to the magnitude of the deposition rate. Growths have been achieved at 
very low temperatures on appropriate substrates if, and only if, the surface 
was relatively smooth and defect free as determined by Kikuchi lines present 
in the RHEED patterns. If the defect density is too high, then epitaxy is 
inhibited by the dominance of growth from the defects. If the deposition 
rate is too high then even with low defect density the growth around a 
defect is so rapid that epitaxy is also limited. Therefore, the understand- 
ing of the growth rate about different types of defects would be helpful in 
assessing the probability of epitaxy for a given system, and we will pursue 
these questions at a later time. 
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In experiment 4, we varied the deposition rate (1 nm/sec **1.8 (IQ 15 ) 
adatoras/cm 2 sec), and figure 13 illustrates the effect of this variation on 
the surface roughness. More experimentation in this area is necessary as we 
will see from the next computer experiment, which also includes variable 
deposition rates. 

In experiment 5, the substrate temperature was varied. Figures 14-17 
illustrate the effect of this variation. The uniformity of the surface is 
progressively improved as the temperature increases and the deposition rate 
decreases, Figure 13 illustrates the roughness factor as a function of time 
and demonstrates the changing surface uniformity. The oscillatory nature of 
the curves is a consequence of the growing multilayers and that surface 
migration tends to fill in the vacancies, ledges and kink sites. The sur- 
face adatoms becomes more mobile at the higher temperatures and tend to fill 
in these sites. This behavior has been observed by Weeks and Gilmer (ref. 

6) for crystal growth from the melt. Figures 14(a-c) illustrate the growth 
sequence at T 3 300 K as compared to figures 14(d-f) which show the growth 
sequence at T a 400 K. Figure 15 illustrates the growth sequence for Rj n 
0.2778 nm/sec and T 53 500 K. Figure 16 illustrates growth occurring by an 
advancing ledge that was formed upon the coalescence of two large clusters. 
Growth by this mechanism has been discussed by Weeks and Gilmer (ref. 6), 
but in this particular case it is statistical in nature rather than the 
dominant phenomena. 

Coalescing clusters at submonolayer levi 8 have been illustrated in 
figure 15 to show the possible development of a grain boundary. Although a 
graphic representation of different growing grain orientations is not easily 

I 

done with the SOJ model, the potential field surrounding a particular defect 
or nucleation site does provide some information on the probable growth 
orientation and, therefore, may permit a way of deciding if the coalescing 
grains will be in, registry, near registry (low angle grain boundary) or 
whether a high angle grain boundary will be formed, 

A plot of nucleation density n s and maximum cluster size N s as 
a function of time is given in figure 19 for a deposition rate of 0.2778 
nm/sec. Clusters of size n£ (i > 1) have a total density of 


y n. = n . (16) 

, l s 
i=z 
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Definite maxima are observed for all temperatures tested. As is apparent 
from the time frames of figures 14-17, the decay in n 8 is due to the 
growth coalescence of clusters. This behavior has also been observed by 
Donahoe and Robbins (Au/NaCl) (ref. 22), Hamilton and Logel (Ag/C,Pd/C) 

(ref. 23), and Corbett and Boswell (Ag/MoSg) (ref. 24). The maximum cluster 
size is also shown to decrease with increasing temperature as w&a also 
observed by Poppa for (Bi/c,Ag/C) (ref. 25). The most probable size of 
clusters for the maxima at T - 300 K and t " 5.5 sec. is approximately 2-3 
adatoms. The initial slope of these curves is the nucle&tion rate which is 
given by 


T 


zo 


* R 


i*+l | E* + (i*+l) AHads - 

KT 


exp 


(17) 


where z is the Zeldovitch factor, o* is the capture number, E* is the 
cluster energy, and i" is the critical cluster size. 

In experiment 6, we studied the effects of annealing. In a similar 
manner to the ordering that occurs in experiment 2, the annealing of a given 
deposition of growth proceeds by surface diffusion or monologue exchange 
between clusters (Ostwald ripening). As discussed previously, cluster 
peripheral motion may also occur to effect an increase in the order of the 
growth. Figure 20(a) illustrates the initial growth condition vised in our 
study. Then with R^ » 0, the substrate temperature is increased. Each 
island or cluster is driven to maximize its number of nearest neighbors 
giving rise to more ordered arrangements as shown in figures for the anneal 
temperatures of T * 600 and t » 700 K. Note that the number of smaller 
nuclei has not noticeably decreased. Figure 19 shows the effect of tempera- 
ture on the average density of the clusters. This same sort of behavior was 

observed by Donohoe and Robins for the Au/NaCl system. Annealing at deposi- 

tion temperature did not seem to have a significant effect even for large 
times. However, when the temperature was increased the low density clusters 
ultimately broke up by nonomer exchange to the larger more stable clusters. 

In the (111) face we assumed that <f> 2 " 0 for ease of computation and 

results for the other crystal faces (111) and (llv/ are not yet available at 

this time. 
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ORIGINAL PAGE IS 
OF POOR QUALITY 

THREE-DIMENSIONAL CONCEPTS AND EXTENSIONS OF MODEL 

/ 

In order to get away from a SOS model one must get involved with the 

three-dimensional geometry of crystal growth. With this purpose in mind, 

consider figure 21 which illustrates a' set of orthogonal axes with basis 
a ^ a a 

vectors 'e‘ 1 =* i, l£> * —• j , “e* 3 * k together with a cubic (P), body 
centered (I) and face centered (F) lattices. 

For the simple cubic lattice all crystal lattice sites are given by 

r = 2£~e'j + 2me^ + 2rie* 3 

where £, m, n are integers. From any lattice point there are six nearest 
neighbor (NN) positions given by the directions ± (n^, n“ 2 , n“ 3 ) where "n 1 = 
2 e i, ng - 2 e~g, lT 3 = 2 e^. There are 12 second nearest neighbors (SNN) given 
by the dxrectxons — (s^, s 2 , s 3 , 84* Sg, 8g) where Sj = n^ + n 2 , 8 g ~ n 2 + 

n 3> s 3 = n l + n 3 > *84 = ”^1 “ n^, ^5 = '^2 “ "n^, Tg = "n^ - n“^. There are 8 
third nearest neighbors (TNN) given by ± (t^, “t 2 , *4) where tTj = iij +n^ 

+ "n 3 , t 2 “"nj - + "n* 3 , £3 = "n^ + ~n 2 "£4 = "n*! “"n 2 "^3* 

For the body centered cubic lattice, all sites are given by 

r = 21 "e^ + 2me^ + 2he^, £,ra,n integers 

together with the set of points 

r - ( 2 l+l)"e^ + ( 2 J+l)e£ + ( 2 K+l)‘e^, I,J,K integers. 

For the neighboring atoms about a given point we have the following direc- 
tions: 

NN directions (8 total) ± (n^, n^, n^, tT^) 
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SNN directions (6 total) ± (s^, s^, 3 * 3 ) 

- 2e^, f 2 - 2e £> s" 3 - 2e^ 

TNN directions (12 total) ± (t’J, T z , T 3 , ’t* 4 , “t 5 , T 6 ) 

tj ■ Sj^ + S £ " 82 ™ Sg 

T 2 - 8^ + r 3 T s - a\ - ^3 

t^-s^ + i^ " si " *2 

Similarly, we find for the face centered cubic crystal that all lattice 
sites are given by 


"r » 2i e*][ + 2me^ + 2ne^, £,m,n integers 

together with the set of points 


"r - (2I+l)eJ + ( 2 J+l)e 2 + 2Ke5 
“r - 21 "ei + (2J+l)e^ + (2K+l)eg 
“r “ (2I+l)e7 + 2J€% + (2K+l)S^ 


I, J,K integers 


For the neighboring adatoms we have: 

NN directions (12 total) ± (n^, rf^, n^‘, iig, ng) 

^ a SJ - e? - ^ 

« 

«2 “ + ^3 *5 " " ®3 

nj - ej + ^*6 " ®2 " ®1 . 

SNN directions (6 total) ± (s^, s^, s*“ 3 ) 

il » 2eJ, aj - 2eJ, s 3 “ 2e^ 



T 

( 

\ 


TNN directions (24 total) ± (t^,' t^, lT 3 , "t^ , . . 



) 


ti « 2e^ + ej + 
t 2 ■ 2e^ + e£ - e^ 

fc 3 “ ®1 + 2S 2 + 
t^ = "e^ + 2e^ + 

t 5 = ej^ + + 2e^ 

*6 = ®l + ®2 ~ 2 ®3 
ty - - 2 e^ + e 2 + e 3 

t 8 “ “ 2 <n + ^2 “ ®3 
tg * ~e*i + 2e£ + 
t 10 = ~e l + 2e£ - 
tn = -e“i + e£ + 2e^ 

*12 = "'®1 + ®2 “ 2e ^ 

For a three-dimensional model of crystal growth we imagine an infinite 
number of lattice sites of one of the three types discussed. We assign some 
initial geometry of occupied sites in the first octant and assume symmetry 
conditions at the boundaries which separate the first octants from the other 
(in order to simplify the upcoming bookkeeping). Each occupied site has a 
potential determined by the number of NN, SNN, and TNN sites which are 
occupied and we write 

U(i,j,k) = NN ijk + SNN.. k * 2 + TNN. jk *3 
Note that the maximum potential when all bonds are in effect are given by: 


SC: U = 6<f>! + 12<f> 2 + 8<f> 3 

BCC: U » 8(f>i + 6(}> 2 + 12t}> 3 

FCC: U = 12<J>i + 6 2 + 24<J> 3 


(18) 
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Here C^) SC does not equal ( ) BCC which does not equal C i ) FCC. 

The symbol <j>^ is used to denote the nearest neighbor bond associated with 
a fixed crystal type. Similar considerations hold for the SNN and TNN bonds 
(J>2, <1)3. 

Having defined an initial geometry we can label those occupied crystal 
sites, where the total bond structure is incomplete, as occupied surface 
sites (OSS). Those lattice sites which are unoccupied, and which are needed 
to complete the bonding of OSS, are labeled unoccupied surface sites (USS). 
Consider only those OSS and USS in the first octant. We can randomly 
deposit adatoms onto USS and convert these sites to OSS or bulk surface 
sites and simultaneously update the potentials at all NN, SNN, and TNN sites 
effected by the deposition as well as recording of the creation of any new 
USS. We can also scan all OSS and assign a random thermal energy to the 
potentials at these sites and determine whether the adatoms at these sites 
remain localized, migrate to NN or SNN sites, incorporate or evaporate. 
Again, the recording of all potential changes of the sites involved, as well 
as the updating of OSS and USS locations, must be performed. 

In order to analyze the three-dimensional growth, the following is 
proposed: Consider the plans illustrated in figure 10 which have direction 

numbers [h]_, h 2 , h 3 ]. This plane passes through the point (la, 0, 0) and 
can be expressed 

hjX + h 2 y + h 3 z = Ih ja (19) 

where I is an fixed integer and h^, h 2 , h 3 are integers. We next con- 
struct the vectors k, ~B*, C which also lie in the plane and are defined as 

A* = -h^ + h 3 e[ 

~B = -h^e^ + h 2 e^ (20) 

C = h 2 e 3 - h 3 e 2 

Let r^ = 2Ie^ be a vector to (la, 0, 0). Then, to determine all lattice 
sites in this plane, we determine rational numbers f, g such that for a 
given crystal structure, 
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r * r Q + fB + gA * xe^ + ye 2 + ze 3 

where x, y, z are lattice sites. In the case of a simple cubic crystal we 
would require that 

2Ie^ + fChje^ - hge^) + ®^1®3 “ ^3^1^ 

where I, hp h 2 , hg, Z, in, n are integers. 

21 “ fh 2 ~ gh 2 ** 2 Z 
fhj = 2m 

ghj “ 2n 

For these equations to possess integer solutions, Z , m, n must satisfy 
hi “ hj£ + h 2 m + h 3 n (21) 

Consider first the case 1=0, then integer solutions ( Z y m, n) to hjA + h 2 m 
+ h 3 n = 0 are given by the points A = (hg, 0, hj), B = (h 2 , hp 0) and C = 
(0, hg, h 2 ). These points define the primitive vectors A, B, C in (20) 
and the primitive cells in figure 22. The primitive cell A, B has a maximum 
of (h^l) interior solutions. Similarly, the primitive cells B, C and A, C 
have a maximum of h 2 ~l and hg-1 solutions. In the primitive A, B cell of 
figure 22, we let n take on the integer values from 1 to (h^-l). For a 
fixed value of n we let m vary from 1 to (h^-1) and find the integer 
solution for the other variable Z. A similar analysis applies to the 
primitive cells B, C and A, C. 

Example : Consider the (295) plane of a simple cubic crystal. Here we want 

integer solutions ( Z , m, n) to the equation 

2Z + 9m + 5n = 0 

We take the obvious solutions (0, 0, 0), A = (5, 0, 2), B = (9, 2, 0) and C 
= (0, 5, 9) and construct primitive cells. For the A, B cell we let n = 1 
and let m = 1 (only one solution) and solve for Z. Thus we develop for 


“ 2£e^‘ + 2me^ + 2ne^ 

This equation requires that 
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° RiGl ^ OR QU AL'TY 


n * T the equation, 2Z + 9m - 5 83 0 with the solution of (7, T, T). Other 
solutions are (7 - 9K, 1 - 2K, T) . K “ ± 0, 1, 2, . . . Similarly, for the 

A, C primitive cell, we let Z ■ 1, 2, 3, 4 and m ■ 1, 2, 3, 4 and we solve 

for n. Thus: 

for Z = 1, 2 + 9m + 5n =» 0 and for m a 1, 2, 3, 4 we find (1, 2, 4) is 

the solution with (1, 2 + 5K, -4 - 9K) as the general solution; 

for Z = 2 (4 + 9m + 5n » 0) + (2, 4, 8) + (2, 4 + 5K, -8 - 9K); 

for Z = 3, 6 + 9m + 5n = 0 + (3 , 1, 3) ( 3 , 1 + 5K, 3 - 9K) ; 

for Z = 4, 8 + 9m + 5n - 0 (4, 3, 7) (4, 3 + 5K, -7 - 9K) ; 

For the B, C primitive cell we let n = I, 2, 3, ?, 5, 6, 7, 8 and Z - 1, 

2, 3, 4, 5, 6, 7, 8 and solve for integer solutions ra. This gives 


n = -1, 
n = -2, 
n = -3, 
n - -4, 
n = -5, 
n = -6, 
n = -7, 
n = -8, 


- 5 = 0 + (7, 1, 1) (7 


2Z + 9m 
2£ + 9m ~ 10 = 0 + 
2Z + 9m - 15 - 0 ♦ 
2£+9m-20=0> 
2Z + 9ra - 25 = 0 -*■ 
2£+9m-30=0+ 
2Ji + 9m-35 = 0 + 
2Z + 9ra - 40 = 0 + 


(5, 0, 2) -*■ (5 
(3, 1,1) (3 

(1, 2, 4) + (1 

(8, 1, 3) (8 

( 6 , 2 , 6 ) > (6 

(4, 3, 7) + (4 

(2, 4, "8) (2 


+ 9K, 7 - 2K, T) 
+ 9K, - 2K, 2) 

+ 9K, 1 - 2K, 3) 
+ 9K, 2 - 2K, 7) 
+ 9K, 1 - 2K, 7) 
+ 9K, 2 - 2K, 6) 
+ 9K, 3 - 2K, 7) 
+ 9K, 4 - 2K, 8) 


The example illustrates that we can determine and count the coverage of 
a given region in a given plane. A similar counting procedure can be 
established for all parallel regions in planes parallel to a given plane as 
in figures lO(b-c). We can then plot coverage vs. distance as we move out- 
ward on parallel planes. Eventually, the coverage becomes -zero and we can 
stop our counting procedure. Those directions for which the coverage vs. 
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distance curves are like step functions, can be labeled as preferred direc- 
tions of growth. Also the distance in these directions (number of growth 
planes measured from some reference) will increase with time at a faster 
rate than other nonpreferred growth directions. During the simulation of a 
three-dimensional crystal growth only the lower ordered planes need to be 
considered as the majority of the growth occurs on these planes. 

CONCLUSIONS 

An SOS Monte Carlo computer simulation utilizing a potential scaling 
technique has been developed to model the initial stages of thin film 
growth. The model makes use of ordered statistics to simulate the surface 
activity of a 20 x 20 array with periodic boundary conditions. Adsorption, 
desorption, surface migration, incorporation and substrate atom diffusion to 
the surface are considered. The results of several computer experiments 
show consistency with the expected behavior of thin-film growth. Surface 
migration data taken at different substrate temperatures returned the acti- 
vation energy as determined by an Arrhonius plot. Dispersed adatoms were 
observed to cluster into dimers, trimers, and finally a single cluster in 
its most stable configuration. A point defect was designed by varying the 
interaction potential in the x, y, z directions and illustrated that 
preferred f,fWt> 7 th occurred at such defects. This suggests a technique for 
studying such defects,, 

The experimentation, with varying deposition rate and substrate temper- 
ature, modeled the expected behavior of thin film growth by nucleation, 
cluster growth and then coalescence of clusters. The nucleation rate was 
proportional to the square of the deposition rate in agreement with the 
atomic theory of nucleation. Finally, annealing experiments showed the time 
ordering of clusters with increased temperatures. 

By adjusting the potential interaction changes and by varying the ini- 
tial geometry and type of surface adatoms, various A-B crystal growth 
models can be investigated. The choice of the A on B, A on A, B on B or B 
on A interaction should produce interesting types of thin film growth 
phenomenon, and this area still has to be investigated. Also three-dimen- 
sional extensions will require large computer times and storage, but the 
basic scanning, testing and updating of potentials will be the main concern. 
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APPENDIX; COMPUTER PROGRAM 


The following computer program is representative of that used in the 
studies of this report. It evolved from earlier models and is currently 
undergoing changes in order to better simulate and model AB, AA, BA, and BB 
types of crystal growth. Graphics output is illustrated in figure 23. 
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OKUK& tV-'-l VJ 

OF POOR QUALITY 


PROGRAM CRYSTAL ( INPUT t OUTPUT » TAPES ■INPUT# TAP E6*0UTPUT>SSETj> 

1TAPE7-SSET) 

COMMON/BLK1/H# EPS(30*30)»ITRAC<30»30»20) 

C0MM0N/PLK2/U0jUM<UE»U0S»UMS»UES»DELT 

CQHM0N/BLK3/XLAM»TSNAP»H»T»TMIG»IEVAP»ICREAT 

CQMM0N/BLK4/FAA<9)*D(9)»FAB(9)»FBA(9)*DS(9)<FBB(9) 

COMMON/ BLK5/L AYER 130> 30) 

C0MH0N/BLK6/Nl>N2jLETA>M2»LAST 

COMMON/ BLK7/PHIl>PHI2>PHI3,PHI0>J>PHSl>PHS2f PHS3»PHS0 

COMMQN/BLKB/ISmCH 

C0MM0N/BLK9/A,AS»EVIB>EVIBS,RD 

COMMON/BLK1Q/LL1*LL2*U1»MAXH»TA»TB 

INTEGER H( 30* 30 ) 

NAHELIST/PARAM/ Hj JjT>DELT»MOjNG»MOS*NQS>PHIO»PHSO» 

1RD»TSNAPjISWTCH>A*AS»TSTQPj liO> UM»UE»UQS»UMS»UES»U 1 

C NAMELIST VARIABLES M»20 IS SIZE OF ARRAY, J»1 ( 100 ) SURF ACE> J »2 (lll)SUR 

C J»3 <110)SURFACE»T«TEMP DEG K# DELT»TIME INTERVAL OF SCAN»ND>MQ MIE 

C POTENTIAL OF GROWING MATL>MOS*NQS»MIE POTENTI ALSFOR SUBSTRATE 

M«20 S % T«500> * DELT«1.0E-* t M0«4 $N0«10 

MOS SNOS-IO S PHI0»2.17 SPHS0«2.17 *TREF«300. STST0P«2. 

RD«200 • S TSNAP* ■ 1 SUP— 3.87 $ UPS— 3.87 S UM— 1. 

UMS— 1 • tUE , D» S UESjOjt > f S &»5 >61 S A S ~ 3 ■ 7 

ISWTCH«0 

C ISWTCH«1 GIVES PICTURES* ISWTCH«0 NO PICTURES IF ISWTCH«0 

C SET R1«0 IN JCL AND I F ISUTCH«1*SET Rl«0 IN JCL C AFTER LGO ) 

C INTERACTION WITH SUBSV l* l T E> PHIO«HEAT EVAP-HEAT ADSORPTION* PHSO IS PHIO 

C OF SUBSTRATE SURFACE*ll)*RATE OF ADATOM ^ EPDSITjADATOMS/SECj 

C A* AS* CRYSTAL LATTICE SPAC INGj.TSTOP«STO- TIME*UO»HEAT EVAP*UM« 

C MIGRATION BARR I ER> ASSUMED CONSTANT*UE«QI * EVAP REFERENCE*U1 NOT USED NOW 

C POTENTIAL ENERGIES IN ELECT. VOLTS U EV*23KCAL/MQLE*M»20 FOR GRAPHICS 

CALL PSEUDO 

500 CONTINUE 

DO 12 I !■ 1* 30 

DO 13 JJ«1*30 : __ 

DO 14 KK«1*20 

ITRAC(II* JJ*KK)«Q 

14 CONTINUE 

13 CONTINUE 

12 CONTINUE 

MAXH-1 

DO B 1 1 ■!* 30 

DO 9 J J ■!> 30 

H<II*JJ)«0 

EPS ( 1 1* J J ) ■0. 

LAYER ( 1 1* J J ) «0 

9 CONTINUE 

e CONTINUE 

READ ( 5# PARAM ) 

IF< EOF ( 5 )") 600*601 











600 

WRITE(6*603) 


603 

FORMAT (1X*28HEND OF FILE ENCOUNTERED-HELP 

5 

CALL CALPLT(0. 0*0. 0*999) 

STOP 6666 

601 

CONTINUE 



LLl® - 


M2«M*M 

LL2«2+H 

CALL SETUP<M0S*N0S*PHS1*U0S*UMS*TREF*PHS0 *J*ALPHA2*ALPHA3) 

P HS2«ALPHA2*PHS1 

PHS3«ALPHA3*PHS1 

CALL SETUP(MD*N0»PHI1*UQ * UM> TREF* PHIO^ * ALPHA2* ALPHA3 ) 

PHI2«ALPHA2»PHI1 

PHI3«ALPHA3»PHI1 

C H«SIDE OF SQUARE ARRAY M HAX«60 

TMIG«Q 

XLAH"RD*DELT 

IEVAP«Q 

ICR EAT»Q 

C INITIALIZE SUBSTRATE GEOMETRY AMD POTENTIALS 

DO 11 J J«1»3Q 

DO 10 1 1 ■ i* 30 _______________ 

10 LAYERC II> J J ) « 0 

11 CONTINUE 

CALL INITIAL 

EHEAN»T*(8.62E-5) 

QD-UM-UQ— PHI 0 

WRITE ( 6* 399 ) 

WRITE (6*400 ) J*DELT*RD*T*M 

WRITE ( 6*401 ) H2 * XLAM » EHE AN* TSN AP»TSTOP 

WRITE (6*402) PHI1* PHI2* PHI3* PHI 0* ALPHA2 

WRITE ( 6*403 ) PHS1*PHS2*PHS3* PHSQ* ALPHA3 

WRITE ( 6* 404) UO*UM*UE* A*U1 

WRITE ( 6*405 ) UPS *UMS*UES * AS* QD 

WRITE (6* 407) HQ* NO* NOS* NOS 

WRITE(6*4Q6) FBB*FAA*D*DS 

407 F0RMAT(T5*3HMQ«*I3*T15*3HN0»*I3*T25*4HM0S«*I3*T35* 4HNQS»*I3) 

399 FQRHAT(1X*15HPARAMETERS ARE* ) 

400 F0RHAT(T10*2HJ«*I3*T34*5HDELT«*1PE14.7*T63* 3HR0»* 0PF7. 1* T91* 2HT«» 

1F8.1*T118*2HM»*I5) 

401 FORMAT (T6* 6H M2 «*I6 *T33*6HXLAM ■* E14 . 6*T60* 6HEHE AN«* F8 .4# T9Q* 

16HTSN AP«* F5. 3* T114* 6HTSt QP«* F7. 1_) 

402 F0RHAT(T7*5HPHI1«*F10.4*T34* 5H) AIZ*> F10. 4*T61* 5HPHI3«* Fjp. 4*T88* 

ISHPHIO^jFIO. 4*T113* 7HALPHA2** F10.4 ) 

403 F0RMAT(T7*5HPHS1»*F10.4*T34*5HPHS2»*F10.4*T61*5HPHS3«*F10.4*T88* 

15HPHS0«»F10.4*T112*7HALPHA3«*F10.4 ) 

404 F0RHAT(T9*3HU0«*F10.4*T36*3HUH»*F10.4*T63*3HUE»*F10.4*T91*2HA«* 

1F10.4 *T115*5HHADS«*F7.3 ) 

405 F0RMAT(T8*4HU0S-*Fl0.4*T35*4Hl‘r,S«,FlO.4*T62*4HUE$«*F10.4*T9O* 
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ORIGINAL PASS IS ! 

OF POOR QUALITY * j 

i 

i 

13HAS»jF10.4 »T117*3HQD«»F7«3) j 

406 F0RMAT(T10»4HFBB»«g9(2X>F8.4) » /> T9» 4HF AA« # 9< 2X* F8»4 ) j /» TIP* 2HD«» j 

19(2XjF8.4)WjT9;3HDS»>9<2XjF 8.4)>/// ) j 

WRITE(6,408)FAB»FBA j 

408 FQRMAT<T10»4HFAB->9(2X ,J_* .4 ) » /> TIP# 4HFBA«j 9( 2X» F8 ,4 ) ) ! 

WR ITE { 6j 308 ) 

308' FORMAT ( 1H1 ) i 

CALL STOF IN< 0»TST0P ) j 

GO TO 500 j 

END j 

SUBROUTINE GROW < ITHREE ) 

CQMH0N/BLK5/LAYER(30f 30) j| 

DIMENSION Xl(4)>X2(4)jYl(4)»Y2<4) j 

SF«20./3« | 

XI ( 3 ) »0 > ; 

Yl(3)«0. 

X1(4)«SF 

Y1 ( 4 ) *SF 

X2 ( 3 ) «0« 

V2(3)»0. 

VOf A\.CC 

AblTI-tfr ( 

Y2(4)«SF 

IFC ITHREE »EQ»1) CALL C ALPLT ( 0. 0, 12. Op-3 ) 

FS-3W20. j 

IFCITHREE«EQ>1)G0 TO 1 „ j 

CALL CALPLT(0.0»->4.0j-3) j 

GO TO 4 j 

X CALL CALPLT(0«0*-3..0»--3) 

4 CONTINUE 

CALL GRID(0»0,0.0*FS»FS>20»20) 

C GRID COMPLETE SHADE IN LEVELS 

DO 10 II«3f22 

D O 15 J J «3» 22 

IH»M0D<LAVER(IIjJJ)»4) 

I F ( IN# EQ • 3 ) GO TO 15 

xim«n-3 

YUl)* JJ-2 

XH2)«II-2 

Yl(2)»JJ-2 

X2 ( 1 ) « 1 1-3 

Y2(l )«JJ-3 

X2m«II-2 

Y2(2)*JJ~3 

C SHADE AREA BTWN POINTS 

INT«1»IH+2*(IH-1)»IH 

CALL HAFTQNE_(XX*Ylg2j>X2»Y2j2»INT) 

15 CONTINUE 

10 CONTINUE 

IFCITHREE.LT.3)RETURN 
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ORIGINAL PAG?; 18 
OF POOR QUALITY 


CALL NFRAME 

ITHREE *0 

RETURN 

END 

SUBROUTINE SETUP < MQjNO* APHI* AUO j AUHj AT» APH , JA# ALPHA2# ALPHA3 ) 

PHIQ«APH 

J «J A 

U0« AUO 

UM«AUM 

T«AT 

R1»1«/SQRT (2> ) 

R2«l « / SORT ( 3 . ) 

R3»FL0AT<N0)/FLQAT<M0> 

A21«(R1»*N0)~R3*(R1**M0) 

A31»(R2*»N0)-R3»(R2»*M0) 

Dl«l.-R3 

A2»A21/D1 

A3»A31/D1 ' 

IF( J «EQ«2 )G0 TO 111 _ 

IFtJ.EQ.3)G0 TO 110 

C 100 FACE 

ALPHA2»A2 

ALPHA3»0» 

PHI1«PHI0/(1.»ALPHA2) 

phiwhiim. ; 

GO TO 12 

111 ALPHA2«A3 

AL PH A3 »0 « ; ___ 

PHIl«PHI0/<6.+2.*ALPHA2) 

GO TO 12 

110 AL PHA2»A2 

ALPHA3«A3 

PHI1»PHI0/(2.*(1.+ALPHA2)»4.»ALPHA3) 

12 CONTINUE 

PHI1— PHI1 

APHWHI1 

RETURN 

END __ 

SUBROUTINE FACE( F»DjA*J*PHI1jPHI2»PHI3»PHI0»EVIB«T) 

DIMENSION F(9>jD<9) 

A2«A/SQRT(2. ) 

A3«A/SQRT (2>/3» ) 

C PHI1>PHI2»PHI3 SHOULD BE NEGATIVE 

C J»l>2 OR 3 

IF(J.EQ.B) 60 T O 110 

IF(J»EQ»2) GO TO 111 

100 DO 16 IJ «1# 4 

IK*2*IJ-1 

IL-2+IJ 
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ORIGINAL PAGE IS 

OF POOR QUALITY 


D ( IK ) ■ A 

D ( IL ) ■ A2 

F ( IK L“PHI 2 

is EijD.sP.nii 

F ( 9 ) »PH I0 

GD TQ 10 

111 DO 17 IJ«2 iA 

D(1J)«A2 

P<IJ+4)«A.2 

F ( I J ) «PHI 1 

17 F < I J +4 ) «PHI1 

D ( 1 ) »A3 

D(5)»A3 

FC1WHI2 

F ( 5 ) ■PH I 2 

F,t..9.),?.P-fcLIH _ 

GD TO 15 

110 P P 1 9. 

IK«2»IJ-1 

D ( IK ) »A 3 

18 F ( IK ) «PHI3 

W J A k « A 

U.\ &JLfA 4 

D ( 6 ) «A2 

om«A 

D(8)*A 

F ( 2 ) ■PHI1 

F (6 ) «PHI1 

F(4)«PHI2 

F(8)*PHI2 

F ( 9 ) *PHIQ 

10 CONTINUE 

EVIB»T»(8.62E-5) 

P<9)»0, 

RETURN 

END 

SUBROUTINE INITIAL 

CQMM0N/BIK1/H>EPS(30>30)>ITRAC(30>3Q>20) 

COMMON/BLK2/UO>UM>US>UOS> IMS >UES> PELT 

C0HM0N/BLK3/XLAH>TSNAP>H>T>TMIG>IEVAP>ICREAT 

C0MM0N/BLKA/FAA(9)»D(9)>FAB(9)>FBA(9)>DS<9)>FBB(9? 

COMMON/ BLK 5/LAYER (30 >30) 

C0MH0N/BLK6/N1>N2jLETA>M2>LAST 

COMMON /BLK7/_P HI 1-aP-HI-2> PHI3»PHI0j J> PHSls PHS2> PHS3> PHSO 

C0MM0N/BLK8/ISWTCH 

C0MMQN/BLK9/A> AS> E VIB> EVIBSr RD 

C0MM0N/BLK1Q/LL1>LL2>U1>MAXH>TA>TB 

C H«ARRAY SIZE 

C INITIALIZE SUBSTRATE GEOMETRY AND HEIGHTS H 

INTEGER H ( 30> 30 ) 
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' Vici-i W 

W* »»** » 

pp P007 1 O n AUTY 


DP 312 IK»1»30 

DO 311 I J ■!» 30 

ITRAC ( IK# I J* 1 ) »0 

311 H(IJ»IK)«1 

312 CONTINUE 

CALL FACE(FAA*DS*AS*J* PHS1* PHS2* PHS3# PHSQ* EVI BS*T ) 

CALL FACE(FBB»D»A# J * PHI1» PHI 2* PHI3* PHIO» E VI B* T ) 

C CHECK SIGNS ON POTENTIALS 

DO 20 1 1« l* 9 

FBA( II ) »FAA( II ) 

! FAB (II ) ■FBB ( II ) 

I 20 CONTINUE 

S F8A(9) iFAA (9 ) 4-UQ-UQS 

i FAB(9)«FBB(9)+UPS-U0 

I DO A1 IKrl>3Q 

DO AO T J« 1* 30 _ 

AO EPS ( I J » IK ) ■UPS 

A1 CONTINUE 

IT YPE-0 

, _C ITYPE«0 SUBSTRATE MATERIAL A-TYPE - 

I C ITYPE°1 GROWING MATERIAL B-TYPE 

C FLAT SURFACE FOR NOW 

C THIS IS THE PLACE TO PUT IN GEOMETRY CHANGES INITIALIZING SUBSTRATE 

I C RADIUS OF CURVATURE EXPERIMENT INITIAL I ZED HERE « EXAMPLE s A LEFT EDGE 

: C DO 300 J J ■!» 2 8 

C 300 CALL ADATOM(FAA»bLl»JJ) 

i C THIS IS THE PLACE TO PUT IN A TRAP NEAR CENTER EPS ( 10* 10 )«U0S-TRAP 

C FAA CORRECTED FOR FIRST LAYER ADDITION ONLY. AFTER. FIRST LAYER. BUILD UE 

' WR XTE ( fc* 100 ) 

100 F0RHAT(1X*25HINITIAL SUBSTRATE HEIGHTS »//) 

\ DO 200 IJ«LL1»LL2 

; WRITE (6# 101) (H(IJ*I)*I»LL1*LL2) 

■ 101 ‘ P0RHAT(10X*20I2) 

200 CONTINUE 

DO 201 I J»LL1# LL2 

WRITE (6# 202) ( EPS ( I J* I ) » I«LL1* LL2 ) 

202 F0RHAT(10X*20F6.2) 

201 CONTINUE 

RETURN 


SUBROUTINE DIFF( II# J J> DMIGjIN* JN) 

COHN ON /BLK1/H*EPS(30*30)*ITR AC (30*30*20) 

CDHH0N/BLK2/U0* UH*UE*UOS*UHS*UES»DELT 

C0HH0N/BLK3/XLAH*TSNAP»H*T*THIG*IEVAP,ICREAT 1 

C0HH0N/BLK4/FAA(9)*D(9)»FAB(9)*FBA(9)*DS(9)*FBB(9) 

C0HH0N/BLK5/LAYER(30»30) 

C0HH0N/BLK7/PHU*PHI2* PHI3* PHIO* J* PHS1* PHS2* PHS3* PHSQ 

C0HHDN/BLK8/ISWTCH 

C0MH0N/PLK10/LL1*LL2#U1#HAXH*TA#TB 
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c 


c 

c 



20 


£ 


400 


401 


402 






rr* ft 4 ' **\ 

j. ' i . ’ > 


OF POOR QUALITY 


}. 


i 

i 

i 

} 


DIMENSION IV(8»2)jIS(8»2)>WW(8)»Q(8) 

INTEGER W ( 8 ) 

INTEGER H ( 30» 30 ) 

IC^O 

RANDOM WALK FROM SITE I TO NI 

DMIG*0. 

FIND AVAILABLE SITES (IF ANY) 

IS ( 8 ) IS LIST OF SAVED AVAILABLE SITES* IC IS NUMBER OF SITES 

DO 10 )! JK* 1» 8 

W ( I J K ) *0« 

CONTINUE 

REMOVE ADATOM 

IJK»H( II> J J) 

IF(IJK.LE«Q)ITYPER*Q 

I F ( I JK «LE « 0) GO TO 20 

ITYPER*ITRAC(II»JJj>XJK) 

CONTINUE 

CALL SUBATOM<II*JJjITYPER) 

CALL XYX 

JJ1*JJ+1 

J JO* J J — 1 

v a * » . « 

Ul'UTj 

110*11-1 

I F < JJ1«GT«22> JJ1*3 

IFUJ0.LT.3) JJ0*22 

IFdll.ST. 22)111*3 

IF (II0«LT«3y 110*22 

FIND AVAILABLE MIGRATION SITES s „ 

IQ*H( 1 1; JJ )+l 

Q( 1) *IQ-H( 1 1 ,» J J 1 ) 

IF (Q (1 ) »LE>0)G0 TO 400 

IC*IC+1 

ISdC>l)*II 

IS dC» 2 ) * J J1 

CALL ?BONDdI>JJljWdC)jJ) 

Q(2)*IQ-H(IIliJJ) 

I F ( Q ( 2 ) • LE »0 ) GO TO 401 

IC*IC+1 

ISdC*l)*IIl 



CALL PBONDdIl>JJjWdC)jJ) 

Q(3)*IQ-H( in JJO) . 

IF<Q(3).LE.0)G0 TO 402 

IC*IC»1 

IS(IC>1)*II 

IS dC> 2 ) * J JO 

CA LL PBOND(inJJO»W(IC )>J) 

Q ( 4 ) ■ I Q-H d 1 0* J J ) 

IF(Q(4) .LE.OJGO TO 40/ 
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or:c Kwt'L pr.c-s u 

OF POOR QUALITY 


IC»IC+1 

CALL PBQND(IIO>JJjW(IC)*J) 

IS_( IC»l )»IIO 

IS(IC<2)"JJ ___ 

403 Q ( 5 ) ■IQ-H( II0> J41 ) 

IF(Q(5) » LE »0 ) GO TO 404 

IC-IC + 1 

IS(IC»1)«II 0 

IS ( 1C» 2 ) ■ J J1 

CALL PB0ND(IIQjJJ1jW(IC)»J) 

404 Q ( 6 ) ■ IQ-H ( IllgJJl) 

IF(Q(6)«LE«0)GQ TO 405 

IC«IC+1 

IS(IC#1)-II1 

IS ( 1C» 2 ) ■ J41 

CALL PBQND(II1»JJ1»W(IC)»J) 

405 Q<7) ■IQ-H( III# J JO) 

IF(Q(7) » L E «0 ) GO TO 406 

IC«IC + 1 

IS < IC»1)«II1 

IS ( ICp 2) JO 

CALL PBONDCIIl* JJO»W(IC)*J) 

406 G i 6 ? "IQ-HC IIO> J JOi 

I F ( Q ( 6 ) » LE .0 ) GO TO 407 

IC-IC+1 

IS(IC>1)»II0 

IS { IC> 2 ) ■ J JO 

CALL PBONPdIO# JJ0>W(IC)»J) 

407 CONTINUE 

JL«4+2*(J-l)-3»( J-l)*( J-2) 

C TEST IC 

IF< IC.EQ.0)G0 TO 500 

C WEIGHTED RANDOM WALK 

WT«0 

IF C IC » EQ • 1 ) GO TO 100 

DO 408 IH«1»IC 

WT«WT+W(IH) 

408 MM < I M ) ■ WT 

DO 409 IH«1pIC 

409 MMtIM)«MW( IMl/MT 

R «UR AN ( 0 ) 

IH»1 1 

411 IF(R.LE.WMUM) 

IFtR.LE.WW(IHn g P TO 412 

IH«IM+1 

IF(IM»EQ«IC)MS*IC 

I F ( I M» ECU IC ) GO TO 412 

GO TO 411 

412 IN-ISCMSpI) 
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0RS05NM- PAC-Ti U 
OF POOR QUALITY 


J H-IS < MS * 2 ) 

G O TO 60 

100 IN-IS(1*1) 

JN-IS<1j2) 

GO TO 60 

500 CONTINbc 

C JUMP TO HIGHER LEVEL ONLY IF RANDOM ENERGY GREATER THAN 

C NEAREST NEIGHBOR BONDS - RANDOM WALK TO NN SITES 

R-URAN ( 0) 

MSITE-R»JL+2 

I F C MSITE ♦ LE » J L ) GO TO 567 

C NO MIGRATION 

IN-II 

JN-JJ 

DMIG-0 

570 CALL ADATOM ( INj4N» ITYPER) ! 

RETURN 

567 CONTINUE 

IV(1*1)-II 

IV(1>2)"JJ1 

TUU.1 »«TT 



I Vt 2» 2 ) -J J 0 _ 

IV<3fl)»IIl 

I V( 3j 2 ) ■ J J 

IV(4»1)-II0 

IV ( 4> 2 ) ■ J J 

I V< 5» 1 ) -I II 

IV(5,2)-JJ1 

I V( 6* 1 ) -I 10 

IV<6j2)»JJ0 

IN-IV(MSITEjI) 

JN-I V( MS ITEj 2 ) 

OMIG-D ( MSITE) 

60 CONTINUE .. 

C TEST FOR MIGRATION OVER BOUNDARIES 

IF( IN»LT »LL1 )G0 TO 70 

IF ( IN«GT« LL2 )G0 TO 71 

IF( JN.LT.LLDGO TO 72 

I F ( JN »GT » LL2 ) GO TO 73 

GO TO 80 

70 CONTINUE 

IF(JN.EQ.2)G0 TO 74 

IF ( JN« EQ« 23 E GO TO 75 

IN-22 

GO TO BO 

74 IH^ZZ 

JN-22 

GO TO 80 

75 IN-22 
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ORIGINS- w 
0F POOR QUALITY 


JN*3 

GO TO 80 

71 IF(JN.EQ«2)G0 TO 76 

IF ( JN« EQ. 23) GO TO 77 

IN*3 

GO TO 80 

76 IN*3 

JN*22 

GO TO 80 

77 IN»3 _ 

JN*3 

GO TO BO 

72 IF(IN«EQ.2)G0 TO 74 

IF(IN.EQ.23)G0 TO 76 

JN*22 

GO TO BO 

73 IF(IN«EQ.2)GQ TO 75 

IP ( I N. EQ. 23 ) GO TO 77 

JN*3 

80 CONTINUE 

CALL ADATOM ( IN* JN* ITYPER ) 

DHIG*D ( 2 ) 

RETURN 

|N0 

SUBROUTINE PBOND ( 1 I* J J * NNB * JCODE ) 

COMM ON /BLK1/H*EPS(30*3Q)*I TRAC (30*30*20) 

CQMHQN/BLK5/LAYER(30*30> 

INTEGER Q(6) 

INTEGER H ( 30* 30 ) 

NNB* 1 

J J1*JJ+1 

JJ0*JJ-1 

II1-1I+1 

110*11-1 

IQ*H ( 1 1* J J )♦! 

IF( JCQOE a EQ> 3 ) GO TO 3 

IF( JCODE. EQ.2)G0 TO 2 

Q ( 1 ) * I Q— H ( 1 1 * J J 1 ) 

Q ( 2 ) ■ I Q-H ( 1 1 1 » J J ) 

Q ( 3 ) « IQ-H (II>JJO) 

0(4) * IQ-H ( 1 10» J J ) 

DO 10 IK*1*4 

IP( Q( IK ) . LE.O )NNB*NNB4-1 

10 CONTINUE ; 

IF (NNB .EQ. 3) NNB *1 

RETURN 

2 CONTINUE 

Q ( 1 ) * I Q— H ( 1 1 * J J 1 ) 

0(2) *IQ-H( I II* J J 1 ) 








ouiumL tree is 

OF POOR QUALITY 




Q ( 3 ) ■ I Q— H ( 1 1 1 » J J) 

Q(4) »IQ-H( II* JJO) 

Q<5)«IQ-H(II0,JJ0) 

q<6)«IQ-H(II0*J,M 

DO 20 IK»1> 6 

IF(Q(IK)»LE«0)NNB«NNB+1 1 

20 CONTINUE 

IF(NNB.EQ»7)NN8»1 

RETURN 

3 CONTINUE 

Q(1)«IQ-H(II»JJ1) 

Q ( 2 ) ■ IQ-H( II* J JO ) 

DO 30 IK»1*2 

IF(Q(IK)«LE.Q)NNB»NNB+1 

30 CONTINUE 

IF(NNB»EQ.«3)-NNB.".l 

RETURN 

END 

SUBROUTINE ADATOM ( I*J*ITYPE) 

CQMMON/BLK1/H*SPS<30*30)*ITRAC(30>30>20) 

CQMMQN/BLK4/FAA(9)*D(9)*FAB<9)*FBA<9)*DS<9)*FBBC9) 

CQMM0N/BIK5/LAYER( 30>30) 

DIMENSION W ( 9 ) 

INTEGER H ( 30* 30) 

IFtH(I*J).LE.0)ISP0T«0 

I F ( H ( I * J ) « LE »0 ) GO TO 80 

C GO AHEAD AND ADD NEW ADATOMS THEN CORRECT BOUNDARIES IF NECESSARY 

IS PQ I«I JR A CU f J j JMlj.il 1 . 

80 CONTINUE 

IF(<ISPOT.EQ.O)»AND.(lTYPE.EQ«Q))GO TO 10 

IF((ISPOT.EQ.O).AND>(ITYPE.EQ.l) )60 TO 20 

IFUISPOT.EQ.l) .AND.(ITYPE,EQ.Q) )G0 TO 30 

DO 40 LM«1*9 

40 W ( LM ) ■FBB ( LM ) 

GO TO 100 

10 DO 50 LM»1*9 

50 W(LM)«FAA(LM) 

GO TO 100 

20 DO 60 LM«1*9 

60 W(LH)-FBAtLM) 

GO TO 100 

30 DO 70 LM*1» 9 

70 W(LM)«FAB<LM) 

100 CONTINUE 



EPS(I* J)-EPS(I*J)+W<99 
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ORIGIN , v/™ m 
OF POC., QUALITY 


C 


2 


£ 


80 


40 


10 

50 

20 ~ 

60 

3 CT 

70 

ICO 


EPS(I“1*J+1)«EPS(I-1»J+1)+W<1) 

EPS<I»J+1)«EPS<I»J+1)+W(2) 

EPS<I + l»J + l)«EPS(I+l#J+l)+Wm 

EPS(I + 1»J)«EPS(I + 1»J) + U(4) 

EPS<I+1»J-1)«EPS(I+1»J-1)+W<5) 

EPSCIj J- l)»EPStI# J-l>+W(6) 

EPS ( I-l» J-l ) »EPS < I-l» J-l ) ♦ W( 7 ) 

EPS(I-1*J)«EPS(I-1*J)+W(8) . 

H(I»J)«H(IjJ)+l 

LAYER(I»J)«LAYER(I»J)+1 

ITRAC(I»J»H(IjJ) )«ITYPE 

IS ADATOM AT AN EDGE? 

TEST1«(I-4)»(21-I) 

TEST2«(J-4)* ( 21-J ) 

IF( ( TEST1 . LT.0).0R.(TEST2.LT.0) )G0 TO 2 

.I F 1 ( TEST1.EQ.0). OR . (TEST2.EQ , 0 )) C ALL EDD(I#J) 

RETURN 

CONTINUE 

CALL EDGE ( Wj I» J ) 

RETURN 

END 


m rnmiM - Iai m i oji lia 4> itype ) 


rnMMnu/ai u i / u.cot / an . am . ttd ap nn.in.4At 

^ w ■ ■ . i w i ■» r H/I.I W ' ^ V7 ^ V * W M. I nwvi t Jur t,v / 


C0HHQN/BLK4/FAA(9)»D(9)>FAB(9)*FBA(9)*DS(9)»FBB(9) 

CDMMQN/BLK5/LAYER(30*30) 

DIMENSION W(9) 

INTEGER H ( 30» 30 ) 

£J3_A.tiE.AP_-A N.D.._S.U B.T R A C T_Q LD ADATONS THEN CORRECT BOUNDARIES 

IK»H< I» J )-l 

IF(IK.LE.0)I0N»0 

I F ( I K. LE .0 ) GO TO 80 

ION»ITRAClI*J>IK) 

CONTINUE 


IF((IQN.EQ.O).AND.( ITYPE . EQ. 1 ) ) GO TO 10 
IF( ( ION.ECUO) .AND. (ITYPE.EQ.O) )GQ TO 20 
IF(<IQN.EQ*l).AND.(ITYPE.EQ.l) ) GO TO 30 

DO 40 LH«1j9 

W(LM)«FAA(L M) 

GO TO 100 

DO 50 LH»1j9 

W(LM)«FBA(LM) 

GO TO 100 

DO 60 LH«1#9 

W(LH)«FAA(LM) 

GO TO 100 

DO 70 LM«1»9 

W ( LM ) »F AB ( LM ) 

CONTINUE 
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ORIGIN ■ " ^ 

OF poo.,' ^^Ay-irv: 


EPS(I* J)«EPS<I,J )— W < 9 ) 

EPSci-ijj+n^EPSu-ijj+n-wm 

EPS<I»J+1)«EPS<I»J+1)-W(2) 

EPS<I+ltJ+l)«EPS(I+l»J+l)-W<3) 

EPS<I+1jJ)«EPS(I+1»J)-WC4) 

EPS ( I+l# J-l) «=EPS( I»l> J-D-Wf 5) 

EPS ( I# J— 1 )«EPS(I*J— D— W(6) 

EPS<I-1*J-1)»EPS<I-1»J-1)-W(7) 

EPS(I-1»J)*EPS(I-1*J)-WS8) „ 

ITRACU*J»H(I# J ) )«-l _ 

H(I»J)»HC I>J)-1 _ 

LAYER<IfJ)»lAYER(I*J)-l 

C IS ADATOM AT AN EDGE? 

TEST1* ( 1-4 )♦ ( 21-1 ) 

TEST2« ( )♦ ( 21-J ) 

IF< (TESTl.LT^l.QR. (TEST2.LT. ODGD TD 2 

1F< (TESri.FO.O).OR, (TEST2.EQ.0) >CALL EDD(I*J) 

RETURN 

2 CONTINUE __ 

CALL ESSE ( Ig J ) 

RETURN 

END 

SUBROUTINE RANDOM (XQjX1»RN) 

C XO IS SEED SUPPLIED BY USER 

INTEGER XO,AA#MP»Xl 

AA«3125 

MP»34359738337 

Xl«MOD(XQ»AA#MP) 

RN«FLOAT(Xl)/FLOAT(MP) 

XO«Xl 

RETURN 

END 

SUBROUTINE EDGECW>I»J) 

COMMON/ BLK1/H, EPS (30, 30), I TRAC (30, 30,20) 

DIMENSION W<9) 

INTEGER H ( 30> 30) 

C ADJUST FDR EDGE EFFECTS CASES WHERE I«3# 22, J«3j 22 

IF ( I , EQ ,3 ) GO TO 20 

IF ( I , EQ, 22 ) GO TO 25 

I F ( J ,EQ,3 ) GO TO 30 

IF( J.EQ.22 )G0 TO 35 

RETURN 

20 CONTINUE 

C TOP EDGE 1-3 

IF( J,EQ.3 ) 60 TO 21 

IF ( J , EQ, 22 1 »0 TO 22 
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ORlGir^i 
OF POO.; Q ! JA 


to 


; EPS(22»J— 1)«EPS(2»J-1) 

EPS(22#J)«EPS(2»J) 

EPS(22# J+1)«EPS(2jJ+1) 

EPS < 23* J-l )«EPS ( 3* J-l l 

EPS(23> J) »EPS(3» J ) 

EP SJ 23»J + :i?*tPS(3»J+l) 

K %2ii> J) «H(3iJ) 

RETURN 

21 EPS(22t22)»EPS(2f 2) 

E PS(2»22)«EPS(2>2) 

EPS(22»2)«EPS<2»2) 

EPSi2#23)«EPS(2»3) 

EPS(22>23)«EPS<2»3) 

EPS ( 22 j 3)«£PS(2»3) 

EPS(22j»>«EPS(2»»> 

EPS<23#A)«EPS(3,A) 

EPS(3>23)"EPS(3»3) 

EPS(23»23)«EPS(3»3) 

EPS<23*3)«EPS(3j3) 

EPS(4>23)«EPS<4»3) 

EPS(4»22)«EPS<4»2) 

EPS<3#22)«EPS(3j2) 

EPSC23>2)»EPS(3>2) 

EPS(23»22j»EP3(3*2j 

H(3#23)«H(3>3) 

Ht23,23>«H(3»3) 

H( 23j 3 ) «H{ 3» 3 ) 

RETURN 

22 CONTINUE 

H(23*22?»H<3>22) 

Ht 23* 2 ) «H t 3» 22 ) 

H( 3» 2 ) »H( 3>22 ) 

EPS(22>21)-EPS(2»21) 

EPS(23,21)-EPS(3,21) 

EPS(2#2)«EPS<2#22) 

EPS(22j2)«EPS(2»22) 

EPSi22#22)«EPS(2»22) 

E?S(22>23)«EPS(2>23) 

EPS(22*3)-EPS(2>23) 

EPS(2<3)«EPS(2*23) 

EPS(23»23)«EPS<3>23) 

EPS<23>3)«EPS<3>_23) _ 

EPS<3>3)«EPS<3»23) 

EPS(23>22) «EPS (3>22) 

EPS<23>2)«EPS(3#22) 

EPS(3»2)«EPS(3»22) 

EPS(4>3)»EPS(4j23) 

EPS(4>2)«EPS(4»22) 

RETURN 
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BOTTOM EDGE 1-22 


I F ( J • EQ • 3 ) GO TO 26 


EP$<2# J-1)-EPS(22# J-l) 


EPS(3#J-1)«EPS<23# 


EPSC3# J)-EPS<23# J 


EPS(2#J)«EPS( 22# J ) 


EPS(3»J+1)«EPS(23#J+1) 


EPS<2# J+1)-EPS(22# J +1 


H( 2# J ) *H( 22# J ) 


RETURN 


CONTINUE 


H(22#23)*H(22#3> 


H(2#3)-H(22#3) 


EPS(21#22)»EPS(21#2 


EPS ( 21# 23 ) *EP 


EPSJ2#4> -EPS (22*4) 


EPS(3#4)«EPS(23#4) 


EPS(3#3)-EPS(23#3) 


EPS(3#23)-EPS(23»3 



















EPS{2#21)«EPS<22»21) 

EPS(22#2)«EPS<22#22) 

EPS(2#2)»EPS<22#22) 

EPS(2#22)*EPS(22»22) 

RETURN 

30 CONTINUE 

IF(I«EQ.3)G0 TO 21 

IF ( I «EQ «22 )G0 TO 26 

H(I#23)»H(I>3) 

EPS < 1-1*23) ■EPS(I-1*3) 

EPS(I*23)«EPS(I*3) 

EPS(I+1*23)»EPS<I+1»3> 

EPS<I-1*22)-EPS(I-1*2) 

EPS(I*22)« a EPStI*2) 

I EPS<I»1*22)«EPS(I+1»2) 

RETURN 

35 CONTINUE 

IF ( I »EQ »3 ) GO TO 22 

IF(I »EQ.22)GD TO 27 

H(I*2)«H(I>22) 

EPS(I-1>3)«EPS<I-1* 23) 

EPS(I*3)«EPS(I»23) 

EPS(I + l;3)“EPS(I+l-23) 

ipsii-I* 2 i«ipsii-:Ljl 2 ) 

EPS(I*2)«EPS(I»22) 

EPS(I+1»2)«EPS(I+1»22) 

RETURN 

END 

SUBROUTINE EPD(IjJ) 

COMMON/ BIKl/H* EPS < 30*30) *ITRAC< 30* 30 >20) 

, INTEGER H ( 30# 30 ) 

IF<I.EQ.4)G0 TO 20 

IF(I.EQ.21)G0 TO 25 

IF ( J «EQ«4 )GQ TO 30 

IF( J »EQ»21)G0 TO 35 

1 20 CONTINUE 

I F ( J « E Q«4 ) GO TO 21 

I F ( J.EQ.2DG0 TO 22 


EPS(23#J+1)-EPS(3# J+l) 


EP$(23#J-1)«EPS(3*J-1) 


RETURN 


EPS (3# 23) "EPS (3# 3) 


EPS(4#23)"EPS{4#3) 


EPS(5#23)“EPS(5»3) 


EPS(23#23)*EP$(3#3) _ 


EPS(23*3)«EPS<3*3) 


23 







ORIGIN- ' **\*JZ 
OF POOii QUAUtY 


RETURN 

22 EPS(3j2) w EPS<3j22) 

EPSt4>2)«EPS<4,22) ^ ' 

EPS<5j2>»EPS<5,22> 

EPS(23,2)»EPS<3»22) 

Emg3/.m«Em3f20) 

EPS <23,2l)-EP$(3,21) „ 

EPS<23»22>«EPS(3#22) 

RETURN 

j2u5 C.Q.MJ.I HUE..... 

IF(J.EQ«4)GO TO 26 

I F ( J > EQ >21 ) GO TO 27 

) 

EPS(2»J)«EPS(22»J) 

EPS(2>J+1)«EPS(22>J»1) __ 

- R .ETU . BM 

26 E PS < 2» 3 ) ■ E PS 1 22 j 3 ) 

EPS<2>4)»EPS<22>4) 

EPS(2#5)«EPS(22,5) 

EPS(2»23)«EPS<22»3) 

E PS ( 20 j 23 ) ■EP$(2Qg3) 

E&SC.21#.23)»EPS(21»3) 

EPS<22j23)«EPS<22#3) 

RETURN 

27 EPS<22»2)«EPS(22>22) 

EPS<21j2)«EPS<21»22) 

EP S(20»2)-EPS<20»22) 

EmgjJ21s£l>S.lE2f.22) 

EPS <2>20E«EPS<22»20) 

EPS(2»21)«EPS(22»21) 

E P SJ-2y.A2 1 » E f.S122f.E2.) 

RETU RN 

30 IF ( I »EQ >4 l GO TO 21 

I F ( 1_»EQ>21 )GQ TQ 26 

EPS(I-1j23) m EPS(I-1j3) 

EPS(I»23)*EPS(I»3> L 

EPS(I+1»23)»EPS<H-1j33 

RETURN 

35 I F ( I > EQ >4 ) GO TO 22 

JFjLI J i.Q ? ..aJ..g-P^.T.9.-.2.7 

£PS(I~l>2)«EPSn-l»22) 

EPS(I*2)«£PS(Ij22) 

EPS(I+1»2)»EPS{ 1+1*22) 

RETURN 

END 

SUBROUTINE STOFIN(TSTART*TSTOP) 

: 0HHQN/BLK1/H»EPS(30j30) jITRACOO* 30,20) 

COMMON /BLK2/UQ* UM*UE*UOS* UMS*UES*DELT 

C0MM0N/Bl.K3/XLAM*TSNAP>N*T*TMIGj»IEVAP*ICREAT 
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ORIGINAL PAGE IS’ 

OF POOR QUALITY 


C0MM0N/BLK4/FAA(9)»D(9)jFAB(9)»FBA(9)jDS(9)»FBB(9) 

COMMON/ BLK5/L AYER ( 30# 30) 

CQMMQN/BLK6/N1»N2»LETAjM2»LAST 

COMMON/ BLK7/-PH1.1#.. PH I2a PHI3# PHIQg 4>PHS1> PHS2»PHS3* PHSO 

C0MMQN/BLK8/ISWTCH 

CQMHQN/BLK9/A»AS>EVIBjEVIBS»RD 

CQMM0N/BLK1C / LL1* LL2» U1j MAXH> TA» TB 

C VALUES OF PARAMETERS COME VIA COMMON STATEMENTS 

INTEGER H ( 30# 30 ) 

DIM ENSION I YMIG(400»3 ) 

DIMENSION IRS IT ( 300# 3 ) ~ 

DIMENSION IRC C 20# 2 ) t IV(401) 

DIMENSION DC C 400 ) » HTS ( 10 ) 

ITIME*1«0 

ITHRE E*! 

IF(ISt‘TCH.EQ.O)GO TO 6 

t ALL GROW ( ITHREE ) 

6 CONTINUE 1 

ITHREE«ITHREE+1 

IRPLUS'-'O 

IMIGk«0 

uo 5u4 

DO 503 II«1^00 

IYMI G( 1 1 ^ J ) *0 

503 CONTINUE 

504 CONTINUE 

TTIME*TST ART I 

UEVAP«UE 

UMIGR-UM 

ICOUNT »0 

TOUT »0 

PRQBC»Q» 

PROBE-O. ' 

PROBMiQ* 

C I BULK IS BULK DIFFUSION COUNTER-REMOVE ADATOMS FROM SURFACE 

I BULK*0 

WTX«2, 

C WT1 IS WEIGHT FOR BULK DIFFUSION 

C QB IS BULK DIFFUSION ENERGY BARRIER 

XKT«T»(8.62E-5) . 

T OUT -TSNAP 

QD»UM-UQ-PHIO 

TAUP»(1»0E-12)*EXP(QD/XKT) 

C TAUD IS THE MEAN TIME BETWEEN HOPS 

C NHOP IS THE NO, OF HOPS IN PELT SECONFAA 

2 XN HOP-PELT /TAUD 

XNSAMP«q.0E12)»TAUD 

XN0FP»(1»0E12)»DELT 

XNSCAN TOUT/DELT 
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QB"WT1*QD 

UBULK-UO+PHIO+QB 


IF(UBULK.GE.UE)UBULK»UE~.001 



X NT OUT ■TOUT-DELT/2 • 


INITIALIZE TIME COU M TER FOR CASE WHERE IJK»Q 

TINIT-DELT/2. ~ 


XXLAM-XLAH/DELT 


XXLAM IS REALLY RD-DEPOSITION RATE 


XIJK-XLAM 


IJK-XIJK 


IJK IS NUMBER OF RANDOM DEPOSITIONS DURING DELT 


IF PELT IS VERY SMALL XIJK IS ALWAYS LESS THAN ONE AND HENCE IJK«Q 

IN CASE I JK"Q GO TO 700 SCAN SURFACE AND PERFORM MIGRATIONS/EVAPORATIONS 
DURING THE DELT INTERVAL AND THEN RETURN TO 1 ' 


UPDATE TINIT IF IJK -0 


F ( I JK «EQ • 0) TINIT-TINIT+DELT 


HAS ENOUGH TIME ELAPSED FOR THE ADDITION OF ONE ADATOM? 


XXIJK»XXLAM* TINIT 


IXIJK.XXIJK 


IF ( I XI JK ,GT.O)IJK«IXIJK 


K.GT.O 


IF DEPOSITION RATE IS ZERO GO TO 700 SCAN SURFACE 


I F ( XL AM .LE. 0.)G0 TO 700 ___ 


I F ( I JK .EQ. 0) GO TO 700 


IRPLUS-IRPLUS+IJK 





ran 


ADD THERMALLY ACCOMODAT 


IRS*IRSIT ( 


JCS-IRSIT(IL,2 


ITYPE-IRSIT(IL,3> 


CALL ADATOM(IRS,JCS,ITYPE> 


CONTINUE 


CONTINUE 


ADATOMS 




CREATION, MIGRATION, EVAPORATION 


NOW DO ALL MIGRATIONS FOR SITES WITH FLAGS 


NSCAN-0 


DO 502 JJ1-LL1,LL2 


DO 501 II1-LL1#LL2 


I UIU1I II ITH tm 


IF(LAYER( III, JJ1) ,LE.O)GO TO 501 


RAN-URAN(O) 


ABC"-ALOG(RAN)/XNOFP 
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ICREAT«ICREAT*IX 


IEVAP-IEVAP+IZ-IB 


IMIGR-IMIGR+IY 


IBULK*»IBULK + IB 






CALL DIFF(IIlgJJljDHIG,»NIl»NJl) 

THIG«THIG+DHIG 

GO TO 505 


DNTINUE 


ITYR-0 


MCKMil 


ITYR»ITRAC(II1*JJ1»H(II1*JJ15) 


CONTINUE 


CALL SUBAT0M(II1»JJ1jITYR) 
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505 CONTINUE 

506 CONTINUE 

TTIHE-TTIME+DELT 

C OUTPUT EVE RY XQUT _S1C_0ND_S 

IF (TTIHE >6E » XNTOUT) GO TO 6 52 

GO TO 1 . 

652 CONTINUE 

XNTOUWTIME+TOUT-DELT/2. 

C BEFORE STOPPING DO AVERAGES FOR THIS PELT RUN 

C RE IN UNITS OF APATOWS/SEC TO CONVERT TO NH/SEC MULTI PL Y .BYl *_1* A1VAQQ. 

C RE«EVAPQRATION RATE 

RE»FL0AT( IEVAP ) /TTIHE 

C RG-AVG GROWTH RATE IN NH/SEC 

RG*A*(FLOAT( IRPLUS--IEVAP-IBULK) ) /FLOAT (M2) 

RG« » 1*RG/TTIHE 

C - RG IN UNITS OF NH/SEC 

C RF»ROUGHNESS FACTOR I EBAR«AVG POTENTIAL WELL 

EBAR-0 

DO 100 II ■LL1» LL2 

DO 101 JJ«LL1*LL2 

EBak-EBARtEPS{II*JJ) .... 

101 CONTINUE 

100 CONTINUE 

E BAR «E BAR/ FLOAT ( M2 ) 

PF»ABS< (UP— EBAR I /U0)*100 

C PF IS ROUGHNESS FACTOR BY WAY OF POTENTIAL DEVIATION FROM NORM 

XN«ICOUNT 

CALL. UPDATE (RF; THE TA*TH2> DC *HTS*J> 

IF(ISWTCH.EQ.Q)GO TO 5 ; 

CALL GROW ( ITHREE ) 

5 CONTINUE 

ITHREE»ITHRSE+1 

N-XN 

ICC»0 

DO 50 II* 1* AO O 

IF(DC( II) «EQ«0)G0 TO 50 

IF ( ICC •GT«20 )G0 TO 50 

ICC*ICC+1 

IKC(ICC>1)*II 

IKC (ICC>2 ) ■DC_M_I ) 

50 CONTINUE 

WRITE < 6*298 )TTI HE* < ( IKC (I I* JJ) * J J*l> 2 ) » I I *1* ICC ) 

298 FORMAT (IX* 6HTTI ME*> F8» 3* 5X*37HDENSITY OF CLUSTERS (SIZE* FREQUENCYJ 

1 »;»T20j(10(lHC>I3»lH»*I3tlH)»2X))W»T20>UQCim»I3*lH»*I3»lH>»2X 

2)) ) 

WRITE (6*299) C L AYER < 3* J J ) * J J«3* 22 ) * ( LAYER < 3*J jl*44«3* 22 ) *N*-^J(NHilP 

299 FORMAT ( IPX* API 1*T60* 6H N* >19 *TB0*5H J ** I2*T100* 5HNH0P** 

1E14.7) 

WRITE (6* 300) (LAYER(4*JJ)*JJ S 3*22)*(LAYER(4*'JJ)*JJ"3*2*;* IHXGR* 
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1ICREAT*IEVAP*IBULK 

300 F0RHAT<10X*40I1>T60*6HIMIGR»*I9*T80*7HICREAT«*I9*T100*6HIEVAP «*I9* 

1T120*6HIBULK»*I6) 

WRITE <6* 301 M LAYER ( 5» J J ) * J J«3* 22 ) * < LAYER (5*«U h< ifJ«3* 22 )* THETA* TH2* 

ITAUDfTA 1 

301 F0RHAT(10X*40Il*T60*9HTHETAm«»F8.4*T8Q*9HTHE.A<2E-*F8.3*T100» 

15HTAUD»*E14.5*T122*3HTA»*F5.0) 

WRITE ( 6* 302 X LAYER (6* J J U J J«3*22 )*{ LAYER ( 6* JJ)* JJ«3* 22)* IRPLUS* 

1XNH0P* TMIG*TB 

302 F0RHAT(10X*40Il*T60*7HIRPLUS»*r9*T80*5HNH0P«*E14»7*T100*5HTMIG-* 

1E14,7*T122*3HTB«*F5.0) 

WRITE (6*303) ( LAYER (7* J J )* JJ-3*22)* (LAYER (7* JJ)* JJ«3* 22)* RE*RG*RF 

303 FORMAT ( 1CX*4QI1*T60* 3HRE«*F9»4*T80*3HRG»*F9« 4*T100* 3HRF»*F9»4) 

WRITE (6*304H LAYER ( 8* JJ)*JJ«3* 22 )* (LAYER ( 8* J J )* J J «3* 22 ) * XNQFP* 

IE BAR* XNS AMP 

304 FORMAT (10X*40Il*T60*5HN0FP«*E14«7*T8Q*5HEBAR a *F8>3»2X» 

110HNSAHPPH0P«*E14«7) 

WRITE<6*310HLAYER{9»JJ)*JJ»3*22>*(LAYER<9*JJ1* JJ«3*22)*HTS 

310 FORMAT <10X»40I1*T6Q*4HHTS“*10F4»0 ) 

DO 201 1 I »10* 22 

WRITE(6*305)(LAYER(II*JJ)*JJ m 3*22)*(LAYER(II*JJ)*JJ*3*22) 

305 ruRnAT(lGX*40Il» 

201 CONTINUE 

IQUT«I0UT+1 

ITEST1»H0D(I0UT*2) 

IFtITESTl.EQ»0)G0 TO 60 

WRITE ( 6*306 ) 

306 FORMAT ( /) 

GO TO 61 

60 CONTINUE 

WRITE ( 6* 311 ) 

311 FORMAT C 1H1 ) 

61 CONTINUE 

IF(TTIME.GE.TSTOP)GQ TO 876 

IF((TTIHE+DELT/2»)»LT« TSTOP) GO TO 1 

876 CONTINUE 

IF( ISWTCH> EQ>1 ) CALL NFRAME 

RETURN 

END 

SUBROUTINE XYX 

C0MMQN/BLKl/H*EPS(30*30)*ITRAC<30*30*20> 

CQMMQN/BLK5/LAYER(30*30) 

C UPDATE FRINGES OF PERIODIC PATTERN 

INTEGER H ( 30* 30 ) 

DO 6 1*3*22 

LAYERt 1* X ) -LAYER < 21* I) 

H ( 1* I ) «H{ 21* I ) 

L AYER ( 2* I ) "LAYER ( 22* I? 

H (2* I ) «H( 22* I ) 
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LAYERt 23* I ) "LAYER ( 3* I) 

H(23*I)»H(3*I) 

LAYER ( 24# I ) ■LAYER < 4* I) 

H(24*I)«H(4* I> 

LA YER ( I* 1 ) -LAYER < 1*21) 

H(I*1)»H<I*21) 

LAYERt I* 2 ) «LAYER ( UZ 2) 

H 1 1* 2 > W H t 1*22) 

LAYERt 1*23) "LAYER ( 1*3 ) 

H<I*23)«H(I*3) 

L AYER ( I*24)«LAYER< I* 4) 

H(I*24)«H(I*4) 

6 CONTINUE 

C CORNERS UPDATED 

LAYERt 1* 1) "LAYER (21* 21) 

H 1 1* 1 ) *H t 21* 21 ) 

LAYERt 1* 2 ) ■ LAYER ( 21* 22 ) 

H 1 1* 2 ) "Ht 21* 22 ) 

LAYERt 2* 1 ) "LAYER C 22* 21 ) 

H ( 2* 1 ) ■Ht_22j> 21 ) 

LAYERt 2* 2)«LAYER(22*22) 

Ht 2* 2 ) »H( 22* 22 ) 

L AYER ( 23* 1) ■ LAYER t 3*21) 

H ( 23* 1 ) »H ( 3* 21 ) 

LAYER 1 23# 2 ) -LAYER ( 3*22) 

H ( 23* 2 ) »H ( 3* 2 2} 

LAYERt 24* 1 ) ■LAYER t 4*21) 

H ( 24* 1 ) »H ( 4* 21 ) 

LAYERt 24* 2 ) "LAYER t 4*22) 

H(24*2>«H(4*22) 

LAYERt 1* 23) "LAYER ( 21* 3) 

H ( 1* 23) »H ( 21* 3 ) 

LAYERt 1*24)-LAYER(21*4 ) 

H 1 1* 24 ) «H ( 21*4 ) 

LAYERt 2*23)«LAYER(22* 3) 

H ( 2* 23 ) »H ( 22* 3 ) 

LAYERt 2*24) "LAYER 1 22* 4) 

Ht 2# 24) «H ( 22* 4 ) 

LAYER(23*23)«LAYER( 3* 3) 

H ( 23* 23 ) «H 1 3* 3 ) 

LAYER(23*24)"LAYER( 3* 4) 

H 1 23# 24 ) "H ( 3* 4 ) 

LAYER (24*23) «LAYERt 4* 3) 

H t 24* 23 ) W H 1 4* 3 ) 

LAYERt 24*24) «L AYER ( 4* 4) 

Ht24*24)«H(4*4) 

DO 7 I «3* 22 

EPS(1*I)«EPS(21*I) 

EPS 1 2* I )«EPS 1 22* I ) 


j 


I 
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roit3i£^/*nri(3>4 


EPS(24j>23)-EPS<4*3) 


EPS(24#24>-EPS<4*4> 


RETURN 


END 


SUBROUTINE URNS(L* IRSIT>M) 


DIMENSION IRSITC 300» 3) 


ALPHA DEPENDENT UPON THICKNESS OF B ON 


AO-.l 


BO-,2 


ALPHA*AO*EXP<-BO*MAXH) 


itaisiTv 


DO 10 1-1*1 


Y1-URAN(0 


Y2-URANC0) 


IRSIT(I#li-LLl+Yl*M 


IRSIT(I#2)-LLH-Y2*M 


IRSIT(I»3)-: 


IF(Yl.lT.ALPHA»IRSIT<I*3)-0 


CONTINUE 


RETURN 
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INTEGER B<30»30)»C(30j30) 

INTEGER H ( 30» 30) 

DIMENSION DC (400)»HTS( 10) 

DIMENSION NEXT ( 400*2 ) 

C RF»ROUGHNESS FACTOR » RF«1 FOR FLAT SURFACE 

C THETA- C DVSRAGE* DC "DENSITY OF C LUSTERS ( 1* 2* 3* ...» 10 ) 

TS-0 

TA«0 

DO 80 II«3»22 

DO 81 JJ«3»22 

IJ »H ( I Ig J J ) 

IF(ITRAC(II»JJ»IJ).EQ.0)TA«TA+1 

IF(ITRAC<II«.JJ»IJ).E0«1)TB»TB+1 

81 CONTINUE 

80 CONTINUE 

DO 5 I»l>400 " 

5 DC ( I ) «0« 

CALL XYX 

C THETA "COVERAGE 

DO 10 II«1»10 

10 HTS(II)«0« 

DO 5Q 183-22 

DO 50 J»3*22 

DO 51 IK«1*10 

IF ( LAYER ( I*J)«GE»IK)HTS( IK ) «HTS ( IK ) +1 

51 CONTINUE 

50 CONTINUE 

THETA»HTSm/400 

TH2«HTSt2) MOO 

MAX«1 

DO 100 . II»1j 10 

IF(HTS(II ) »GE.20, )MAX-II 

100 CONTINUE 

MAXH»M AX 

DO 8 1*1*24 

DO 8 J«1 *24 

C(I> J)«l 

B(Ij J)»0. 

8 CONTINUE 

DO 9 I«l*24 

DO 9 J«1j24 

IF<LAYER(I»J).GE.MAX)B<I».M«1. 

9 CONTINUE 

SUH«0 » 

DO 60 I»3»22 ; 

DO 60 J»3»22 

60 SUM«SUM+IA8S(H(I»J )-H( I* J + l) ) 

RF«1+SUMM00. 

ILST-0 
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INXT-0 

ICT»0 

DO 20 I«3»22 

DO 20 J«3>22 

IF(B(I» J) .EQ.0)G0 TO 20 


LJ-J j 

CALL CORRECT(LIjLJ) j 

25 ICT«ICT+1 ] 

ILJ«LJ + 1 j 

IF(B(LI»ILJ).EQ.O)GO TO 21 f 

IF(C(LI,ILJ)«EQ.O)GO TO 21 ’ 

ILST »ILST»1 

NEXTtl t ST# 1)«LI 

NEXT(XlLST»2)«ILJ 

CALL CORRECT ( LI# IL J ) 

21 ILJ«LJ-1 

IF(B(LI,ILJ) »EQ.O)GO TO 22 

IF(C(LI,ILJ).EQ.Q)GQ TO 22 

ILST»ILST+1 

NEXT ( ILST# 1) »LI 

NEXTCILST* 2)»ILJ 

CALL CORRECT ( LI» IL J ) 

22 IF(JCQDE«EQ«3)GQ TO 24 

ILI«LI+1 

IF(B(ILI#LJ)»EQ«0)GQ TO 23 

IFCCCILI^LJ) « EQ «0 ) GO TO 23 

ILST«ILST+1 ; !j 

NEXT(ILSTj1)«ILI 

N EXT ( ILST» 2) «L J 

CALL CORRECT ( ILI> L J ) ! 

23 ILMH 

IF<B<ILXjLJ)«EQ.Q)GO TO 74 

IF(C(ILIjLJ)«EQ.O)GO TO 74 

ILST*ILST+1 

NEXT(ILST#1)«ILI 

NEXT ( ILSTf 2) «LJ 

CALL CORRECT ( ILI>L J ) 

74 IF(JC0DE«EQ«1)G0 TO 24 

ILI»LI-1 

ILJ»LJ-1 

IF(B(XLIjILJ)«EQ.O>GO TO 76 

IF(C(ILI>ILJ).£Q.O)GQ TO 76 1 

ILST«ILST + 1 

NEXT(ILSTj1)«ILI 

NEXT(ILST»2)«ILJ 

CALL CORRECTULIjILJ ) 

76 ILI"LI+1 

ILJ-LJ+1 
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IF(B(ILItILJ).EQ.O)GO TD 24 

IF(C(ILItILJ).EQ.O)GO TO 24 

ILST «ILST +1 

NEXT ( ILST t 1 ) ^ILI 

NEXT(ILSTt2)«ILJ 

CALL CORRECT ( ILIt IL J ) 

24 B (LItL J )»Q 

C (Lit LJ ) »Q 

IF ( L I « GT • 3 ) GO TO 30 

B (23t LJ ) *0 

C C 23* L J ) *0 

30 IF(LI.LT.22)G0 TO 31 

B (2tlJ ) *0 

C ( 2tL J ) «0 

31 IF(LJ.GT.3)G0 TO 32 

B(LIt23 )»0 _ 

C ( L It 23 ) *0 

32 IF(L J >LT»22)G0 TO 33 

B ( Lit 2 ) «0 

C(LIt2)«0 

33 CONTINUE 

I F ? IN XT» E w. IL 5T ) GO TO 26 

INXT-INXT+1 

LI«NEXT ( INXTt 1 ) 

L J«NEXT ( INXTt 2 ) 

IF(LI.LT.3)LI»22 

IF(LI.GT.22)LI«3 

IFCLJ.LT. 3)LJ«22 

IF ( L J . GT. 22 ) L J»3 

GO TO 25 

26 ILST«0 

INXT »0 

DC(ICT)»0C(ICT)»1 

ICT»0 

20 CONTINUE 

RETURN 

END 

SUBROUTINE CORRECTCIIt J J ) 

C0HHQN/BLK5/LAYER(30t30) 

CQMHQN/BLK12/C 

INTEGER C ( 30t 30) 

C(IIt JJ)«0 

IF ( J J. EQ» 2>G0 TO 2 

IF( JJ.EQ.3 )GQ TO 3 

IF( J J.EQ«22)G0 TO 22 

IF ( J J . EQ .23 ) GO TO 23 

I F ( I I . EQ .2 ) GO TO 32 

IF ( I I . EQ . 3 )GQ TO 33 

IF ( I I . EQ. 22 ) GO TO 34 
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IF ( 1 1 » EQ« 23 ) GD TD 35 

RETURN 

2 * C(II»22>"0 

17 IF(II.EQ.2)C(22»JJ)«Q 

IF(II.E0.3>C<23# JJ >«0 

IF<II.EQ«22)C(2*JJ)«0 “ 

IF(II*EQ*23)C(3#JJ)"0 
GO TO 24 

3 C<II>23)«0 

GO TO 17 

22 C ( 1 1» 2 ) »Q 

GO TO 17 

23 C(II>3)«0 

GO TO 17 

32 C (22* JJ )*0 

18 IF(JJ.EQ»2)C(IIj22)»Q 

IF(JJ.ECU3)C(II»23)»0 

IF(JJ.EQ.22)C(IIj2)«0 

IF(JJ»EQ»23)C(IIj3)»Q 

GO TO 24 

33 C(23>JJ)«0 

GO TO 18 

n J. *+§*%. I I t . A 

31 UlCfUUJ-V 

GO TO 18 

35 C (3j J J ) «Q 

GO TO 18 

24 CONTINUE 

IF< (II»£Q.3).AND.UJ«EQ.3) )C(23>23>»0 

IF((II.EQ.3).AND.(JJ.EQ.22))C(23j2)«0 

IF( (II.EQ.22) .AND. ( J J , EQ . 3 ) ) C ( 2> 23 ) ■ 0 

IF( (II .EQ.22) .AND, ( JJ.EQ.22) >C(2j2)» 0 

RETURN 

END 

SUBROUTINE MIXUP(ISjNS) 

DIMENSION L ( 401 ) » IS ( 401 ) 

ioi ; 

IFCNS.EQ.DGO TO 100 

DO 50 KS«1»NS 

50 L ( KS ) ■KS 

M»NS 

500 IF<M.EQ.2)G0 TO 102 

R»RANF(0) 

I»1»R*M 

IF( I.EQ.DGO TO 51 

IF(I»EQ.M)GD TD 52 

IS ( IC ) «L ( I ) 

L ( I ) »L ( M ) 

H»H-1 

IOIC+1 
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ADATOMS 


Figure 1. SOS model for crystal growth. 
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Figure 




Profile of surface potentials (a) uniform homogenous surface 
(b) vacancies, dislocations, kinks, jogs and impurities disturb 
the surface potential. 
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Figure 5. FCC model and potential changes associated vith different crystal 
orientations . 
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Figure 7. Flowchart for the Monte Carlo SOS computer simulation of thin 
film growth. 
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Figure 8. A 20 x 20 square array with periodic extension. 
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Figure 10. Clustering of dispersed adatoms on a uniform surface: 
T - 600 K, (a) t - 0, (b) t - 0.1 s, (c) t - 0.3 s, 
(d) t - 0.5 s, (e) t - 0.7 s, (f) t =» 1.0 s. 
























ORIC Ai r 
0F POu.'i Q'jAL.ry 





Figure 12 Growth around a trap: T * 550 K, (a) t 3 0.5 s, (b) t * 
0.8 3, (c) t 3 1.1 s, (d) t 3 1.4 s (e) t 3 1.7 s, 

(f) t 3 2.0 s, 3 0.2778 nra/sec. 
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Figure 14. Thin film growth for * 0.2778 nm/sec for two different 

temperatures, T • 300 K, (a) t * 0.5 s, (b) t * 3.0 s, 

(c) t * 6 .0 s , T * 400 K, (d) t - 0.5 s, (e) t * 3.0 s 
(f) t - 6.0 s. 
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Figure 16. Thin film growth for ■ 0. 5556 ntn/sec, T * 600 K 

fa) t ■ 0.5 s, (b) t ■ n.S s. (c) t * 1.1 s. (d) t 
s, (e) t * 
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Figure 17. Hiin film growth for Rj * 1.389 nra/sec, T ■ 600 K, 

(a) t * 0.25 s, (b) t * 0.4 s, (c ) t ■ 0.55 s, (d) t * 
0.7 s, (e) t * 0.85 s (f) t - 1.0 3. 
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Figure 22. Crystal planes (a) their orientation and 
(b) and (c), and (d) A,B, C,A and B,i 


coverage with distance 
primitive cells. 
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